Small, electrically driven unmanned aircraft are likely to suffer from inferior endurance compared to their larger counterparts, attributable mostly to the limited aerodynamic efficiency of small wings of moderate aspect ratio and the comparably low energy density of available civil battery technology. Upwash exploitation by tight formation flight, inspired by migratory birds, as well as aerial recharging are the most promising controldriven approaches to mitigate this disadvantage while simultaneously opening up new possibilities such as distributed payloads and mission-level redundancy.

List of Figures

1 Flock of grus americana in close formation flight . . . . . . . . . . . . . .

Flight experiments with manned aircraft as well as wind tunnel data suggest that the region in the wake of a predecessor aircraft where a significant reduction in energy consumption can be achieved is however both laterally and vertically of the magnitude of a mere fraction of a wingspan. The docking procedure involved in exchanging batteries between Unmanned Aircraft Systems (UAS) in flight implies similar, if not even tighter, guidance performance requirements. Precise formation flight has, furthermore, to be performed under the adverse influence of only approximately known wake disturbances.

Continuous time higher order sliding mode control (CTHOSM) has been considered as a candidate for this challenging open problem and was successfully applied to simple kinematic models in simulation, where excellent relative position tracking performance can be demonstrated. In this work we study the implications of the presence of inner loop dynamics and discrete implementation at moderate sampling rates and we find that it precludes the application of CTHOSM control to fixed-wing UAS. We propose a predictive discrete sliding mode guidance scheme to approximate the performance of CTHOSM control assuming realistic fixed-wing UAS dynamics. Comparisons with time sampled continuous-time sliding mode control previously proposed show vastly improved chattering and position tracking performance. We show that the proposed guidance scheme in combination with inner load factor tracking loops and a disturbance observer allows for relative position tracking performance compatible with the requirements of upwash exploitation. We address the scalability of local feedback sliding mode control schemes and point out fundamental limitations of sliding mode control with local feedback information when it comes to string stability. To evaluate the performance of the proposed guidance scheme and to facilitate the testing and comparison of new, improved formation flight control techniques, an openly accessible benchmark problem is proposed in this work.

Finally, the closeness of UAS in tight formation flight scenarios makes collision avoidance a major concern. However, existing probabilistic localization algorithms cannot provide guaranteed confidence regions of the relative position of other members of the formation. As an important building block for safe tight formation flight, we present a set membership filter that provides ellipsoidal regions guaranteed to contain the relative positions of the other UAS. Being based on observations of single-frequency Global Satellite Navigation System (GNSS) receivers and ultra-wide band ranging devices, it is compatible with the hardware constraints of small low-cost UAS. Simulations suggest computational efforts compatible with the computational resources typically available onboard small UAS.

Résumé

Les mini-drones à propulsion électrique sont susceptibles d'avoir une endurance inférieure à celle de drones plus grands, en raison, principalement, de l'efficacité aérodynamique limitée des petites ailes de faible allongement. La capacité de stockage limitée des batteries embarquables sur les mini-drones réduit également l'endurance totale. L'exploitation des interactions aérodynamiques, inspirée par les oiseaux migratoires, ainsi que le ravitaillement en vol , sont des approches prometteuses pour améliorer l'endurance des mini-drones tout en permettant une distribution de la charge utile.

Dans le contexte du vol en formation classique des avions, on a observé qu'une réduction significative de la consommation d'énergie de l'avion suiveur peut être obtenue dès lors qu'il se place dans les tourbillons de sillage du prédécesseur. Dans notre contexte, cela implique des déplacements très précis du mini-drone suiveur dont la position relative (latérale et verticale) doit rester à une fraction d'envergure du mini-drone prédécesseur. Par ailleurs, la procédure d'amarrage utilisée dans l'échange de batteries entre drones en vol implique des exigences de performance de guidage similaires, voire plus strictes. Ces fortes contraintes de performance s'accompagnent en outre d'exigences de robustesse particulièrement élevées sur les lois de commande en raison des perturbations mal connues induites par les turbulences de sillage.

La commande par modes glissants d'ordre supérieur en temps continu (CTHOSM) a été considérée comme un candidat prometteur à ce problème ouvert difficile et a été appliquée avec succès à des modèles cinématiques simples. Des performances excellentes de suivi de position relative ont en effet pu être démontrées. Dans nos travaux, nous étudions les implications de la présence de la dynamique de la boucle interne et de l'implémentation en temps discret à des taux d'échantillonnage modérés et constatons alors que l'application de la commande CTHOSM devient impossible. Nous proposons donc un schéma de guidage prédictif discret par modes glissants pour approximer les performances de la commande CTHOSM pour une dynamique réaliste du drone. La comparaison avec la commande par modes glissants en temps continu échantillonnée proposée dans la littérature montre des performances de suivi de position considérablement améliorées. Nous montrons que la loi de guidage proposée, complétée par des boucles de suivi du facteur de charge interne et d'un estimateur des perturbations permet de garantir une précision sur le suivi de position relative compatible avec les exigences de l'exploitation du flux d'air généré par le prédécesseur. Dans nos travaux, nous abordons aussi, avec la notion de « string stability », la problématique liée à la propagation des erreurs de suivi dans des systèmes de commande par modes glissants sous retour d'états locaux et soulignons les limitations fondamentales. Pour évaluer la performance de la loi de guidage proposée et permettre aussi à la communauté de tester diverses améliorations des techniques de commande de vol dédiées au vol en formation de drones, un problème de référence ouvert a été proposé dans le cadre de ce travail.

Enfin, la proximité des drones dans les scénarios de vol en formation serrée nécessite la mise en place de stratégies efficaces d'évitement de collisions. Cependant, les algorithmes de localisation probabilistes existants ne permettent pas la caractérisation de régions de confiance garanties de la position des autres membres de la formation. Dans ce contexte, nous proposons un nouveau filtre ensembliste caractérisant de telles régions de confiance sous forme ellipsoïdale. La mise en oeuvre du filtre proposé repose uniquement sur des mesures accessibles via des récepteurs GNSS à faible coût et des capteurs de distance relatifs. Nos premières évaluations ont montré que les efforts de calcul induits par cette mise en oeuvre restent parfaitement compatibles avec les contraintes des systèmes avioniques des petits drones.

Introduction Background and motivation

Migratory birds routinely perform long range flights, and impressive nonstop distances of more than 8000 km over a flight time of over nine days are reported [START_REF] Gill | Extreme endurance flights by landbirds crossing the pacific ocean: ecological corridor rather than barrier?[END_REF]. This kind of endurance has so far been elusive for man-made unmanned aircraft of similar size. At the same time, the class of small, low cost unmanned aircraft systems (UAS) enjoys a rising market success as their growing maturity and cost efficiency makes them a viable tool for a range of industries, from crop surveillance in agriculture over short-range cargo delivery to mining exploration. Commercial small electrically driven fixed wing UAS, however, typically are able to stay in the air for not more than a few hours. Research efforts focusing on increasing propulsion efficiency and improving onboard energy storage density have succeeded in pushing this envelope in prototypes to several days of endurance over the last years by employing liquid hydrogen fuel cell technology [START_REF] Swider-Lyons | Liquid hydrogen fuel system for small unmanned air vehicles[END_REF]. Other approaches rely on exploiting exogenous energy sources by adding solar cell arrays [START_REF] Oettershagen | High-Fidelity Solar Power Income Modeling for Solar-Electric UAVs: Development and Flight Test Based Verification[END_REF] or localizing thermal updrafts [START_REF] Stolle | A vision-based flight guidance and navigation system for autonomous cross-country soaring UAVs[END_REF]. Tight formation flight has attracted attention over the last decades as a complementary exogenous energy source inspired by nature. Many bird species use efficient formation patterns that allow them to recover part of the energy stored in the upwash generated by the birds in front of them and achievable energy savings of the order of 10 % are reported in the literature [START_REF] Weimerskirch | Energy saving in flight formation[END_REF]. As early as 1990 the feasibility of this approach has been Figure 1: Flock of grus americana in close formation flight demonstrated in humanly piloted flight experiments with manned aircraft [START_REF] Beukenberg | Aerodynamics, performance and control of airplanes in formation flight[END_REF], showing excellent agreement with theoretical predictions. From the NASA AFF program, about a decade later, maximum fuel savings of 18 % of fuel flow for the follower are reported [START_REF] Jake | F/A-18 Performance Benefits Measured During the Autonomous Formation Flight Project[END_REF] for a pair of F/A-18 fighter aircraft under automatic control, further confirming the feasibility of tight formation flight (TFF) for significant energy gains and thus range enhancements. Subsequent more recent flight experiments [START_REF] Pahle | A preliminary flight investigation of formation flight for drag reduction on the C-17 aircraft[END_REF] with C-17 transport aircraft suggest average fuel flow savings of 7 to 8 %. While potential benefits are thus considerably smaller than those achievable by improved energy storage, tight formation flight (TFF) requires only readily available consumer grade commercial-off-the-shelf (COTS) sensors for precise localization and very limited, if at all, airframe modifications. On top of range benefits, the TFF capability is an enabling technology that can unlock a number of other applications with the potential to enlarge the operating envelope of small electrically driven UAS. Aerial recharging and aerial docking are two examples that are not possible today due to a lack of available technologies for safe precise navigation in the close vicinity of another UAS. Position tracking performance requirements of aerial docking obviously depends on the docking mechanism catching the other UAS, but it seems reasonable to assume that they are of the same order of magnitude or stricter than those for upwash exploitation. Any advancement made in the field of small UAS towards safe tight formation flight techniques are of obvious great interest to commercial aviation and could bring an appreciable leap in fuel efficiency to commercial air transport. Small UAS provide in this context an inexpensive testing environment of low complexity that allows to rapidly evaluate the maturity of new techniques. The benefits of upwash exploitation come with supremely challenging guidance performance requirements. Both theoretical considerations [START_REF] Hummel | Aerodynamic aspects of formation flight in birds[END_REF] and the flight experiments of [START_REF] Jake | F/A-18 Performance Benefits Measured During the Autonomous Formation Flight Project[END_REF] suggest that significant upwash gains can only be expected within a small region spanning roughly about 20 % of the wingspan vertically and laterally. Longitudinal separations are less critical due to slow vortex decay. For small UAS this translates into decimeter-level relative position tracking performance and corresponding accuracy in relative localization between UAS. The resulting small separations between UAS make collision avoidance a major concern. Flight safety has gained even more importance with the ongoing integration of civil UAS into National Airspaces (NAS). Despite recent advances, important elements are missing that we consider essential to make tight formation flight a safe mode of operation. In this work, three principal issues are addressed: position tracking performance, the lack of benchmark problems and guaranteed localization under observation errors.

Open problems and contributions

Open problem I: position tracking performance

Reported experimental results for heterogeneous [START_REF] Park | Avionics and Control System Development for Mid-Air Rendezvous of Two Unmanned Aerial Vehicles[END_REF] and homogeneous 2-formations [START_REF] Gu | Design and flight testing evaluation of formation control laws[END_REF][START_REF] Rice | Control performance analysis for autonomous close formation flight experiments[END_REF] of small UAS feature relative guidance errors of the order of meters. In order to enable meaningful range enhancement for small UAS, significant advances need thus to be made.

The performance gap Considerable interest has appeared in the control community over roughly the last two decades to tackle this challenging control problem. Looking at the landscape of existing research, it becomes apparent that two quite different approaches have been taken for outer loop control (at the inner loop level we find mostly attitude tracking controllers, an exception being the use of Nonlinear Dynamic Inversion (NDI) by the authors of [START_REF] Gu | Design and flight testing evaluation of formation control laws[END_REF]). At one end of the theory-application spectrum, established control techniques such as LQ control [START_REF] Park | A new nonlinear guidance logic for trajectory tracking[END_REF] or in-flight-tuned proportionalintegral-derivative (PID) control [START_REF] Wilson | Vision-aided guidance and navigation for close formation flight[END_REF] have been applied and successfully flight tested (see [START_REF] Gu | Design and flight testing evaluation of formation control laws[END_REF], and [START_REF] Rice | Autonomous close formation flight control with fixed wing and quadrotor test beds[END_REF] for a more recent implementation of the same control laws on an electric fixed wing UAS). As mentioned, the achieved position errors exceed the tolerances necessary for significant wake energy gains.

On the other hand, more advanced techniques such as Higher Order Sliding Mode control (HOSM) have been applied to heavily simplified kinematic representations of the relative position tracking problem in simulation, and very successfully so in terms of relative position errors, see [START_REF] Galzi | Formation flight strategy and control using higher order sliding mode[END_REF][START_REF] Yamasaki | Sliding mode based integrated guidance and autopilot for chasing UAV with the concept of time-scaled dynamic inversion[END_REF][START_REF] Galzi | Unmanned rotorcraft tight formation flight control using sliding mode control driven by sliding mode disturbance observers[END_REF][START_REF] Yamasaki | Separate-channel integrated guidance and autopilot for a path-following UAV via high-order sliding modes[END_REF][START_REF] Galzi | Closed-coupled formation flight control using quasicontinuous high-order sliding-mode[END_REF][START_REF] Rao | Sliding mode control-based autopilots for leaderless consensus of unmanned aerial vehicles[END_REF]. Sliding mode control exhibits excellent disturbance rejection properties in continuous time for a very specific class of input-affine systems with matched disturbances. Real-world aircraft formation control violates both central assumptions -input switching in continuous time and the matching condition -that form the very basis of SMC techniques. One approach (see e.g. [START_REF] Galzi | Formation flight strategy and control using higher order sliding mode[END_REF]) to circumvent this fundamental issue is to simplify the relative tracking dynamics until the matching condition is satisfied, corresponding to some form of double-integrator dynamics, raising however questions about the applicability of the thus defined control to real fixed-wing aircraft featuring bandwidth-limited rotatory dynamics and actuators.

Contribution I: Predictive sliding mode guidance law

It is considered a prerequisite to successful tight formation flight to close the gap between insufficient position holding performance reported from flight experiments and the superior performance HOSM control can demonstrate on heavily simplified approximations of the relative position tracking problem. In chapter 2 we present Pre-dictive Discrete Sliding Mode (PDSM) guidance laws that enable position holding performance compatible with the tight requirements of TFF in the vertical channel on a realistic benchmark. We show that, in simulation, PDSM control can provide performance comparable to super twisting sliding mode control for fast sampling and fast inner loops, and vastly superior performance under more realistic conditions. A key ingredient of the control strategy studied in this work is an inner/outer loop interface based on load factor tracking. We provide a prototype LQ vertical load factor tracking control law employing Direct Lift Control (DLC), see section 2.4.4.

Open problem II: benchmark problems

Binary properties of a system under automatic control such as robust stability can usually be evaluated unambiguously for a given system. In contrast to that, control performance is intimately tied to the problem that is used to evaluate a certain control approach. An important obstacle to the systematic and performance-driven advancement of tight formation flight control as of today is the absence of agreed-upon benchmarks that allow to evaluate competing approaches and to make an informed judgment about performance and implementability. Especially for control techniques involving highfrequency control inputs, differences in the chosen modeling of aircraft dynamics and more subtle details such as the choice of simulation sampling times can make the difference between superior control performance and total failure for the same control law, as is elaborated in more detail in chapter 2.

Contribution II: openly available formation flight benchmark

We therefore propose in chapter 1 a benchmark accompanied by an openly available reference implementation. Making available an implementation is expected to significantly lower the threshold towards adopting the benchmark by other researchers, taking into account the corresponding implementation effort of several man-weeks (according to our experience).

Open problem III: guidance safety guarantees

Trajectory planning in the vicinity of other UAS needs to take into account both localization errors and UAS guidance errors to avoid inter-vehicle collisions as well as entering regions in the wake where roll control saturation can occur. The most popular localization filtering strategies based on the Unscented (UKF) [START_REF] Gross | Robust UAV relative navigation with DGPS, INS, and peer-to-peer radio ranging[END_REF] or Extended Kalman filter (EKF) [START_REF] Williamson | An instrumentation system applied to formation flight[END_REF], provide confidence measures by estimating the state estimate covariances that can be used to approximate performance bounds. Their stochastic nature preclude however that collision avoidance guarantees can be given.

Contribution III: Set membership localization filter and constrained sliding surface

In chapter 3 we describe a set membership localization filter that computes guaranteed bounds on the position of another UAS based on GNSS carrier phase and UWB range observations. In chapter 2 we exploit the guaranteed relative position sets by designing an integral discrete sliding surface using standard Model Predictive Control (MPC) techniques that enables guaranteed collision avoidance in sliding mode.

Chapter 1 A formation flight benchmark

Benchmark problems play an essential role in every performance-driven field of research, as they allow to quantitatively compare different approaches and to overcome the problems of solution-specific modeling 1 . As of today, there is no lack of contributions that hold the promise of enabling the performance leap necessary for tight formation flight, a lack of reference problems is however evident. This section presents a benchmark problem for automatic aircraft formation flight [START_REF] Bolting | An open benchmark for distributed formation flight control of fixed-wing unmanned aircraft systems[END_REF]. It provides a heavily simplified kinematic model including first-order load factor tracking loop approximations as well as a realistic six degrees of freedom nonlinear model of an electric glider UAS including servo actuator and engine models. A 3D benchmark trajectory and wind and turbulence modeling put guidance laws through a sequence of tracking problems of increasing difficulty, from a straight cruise flight section over climbdescent maneuvers through a level turn into successive helical climb-descent-turn maneuvers. To be useful, a benchmark must be accessible for other researchers, and to this end a Matlab ® /Simulink ® implementation has been made accessible as an open-source project (Formation Flight Benchmark (FFB) [START_REF] Bolting | Formation Flying Benchmark[END_REF]).

Standard Vehicle Model

Every controller design model inevitably involves simplifications that make it a more or less rough approximation of the actual process to be controlled. The performance of various control techniques applied to the problem of tight formation flight is intimately related to the set of simplifying assumptions adopted for modeling the system dynamics. As an example, excellent performance can be achieved by second order sliding mode control when assuming double integrator dynamics and high controller and 1 By this we understand choosing a problem that is tailored to a specific solution approach and inherently favors its success. One example is the evaluation of techniques involving high-frequency control inputs such as sliding mode control on plant approximations that lack important dynamics model sampling rates (see e.g. [START_REF] Galzi | UAV formations control using high order sliding modes[END_REF]), i.e. conditions close to the controller design model. Sampling times compatible with avionics hardware constraints and the introduction of inner loop tracking dynamics lead to unacceptable chattering in the input channel and unacceptable position errors, see sc. 2.4.4. To systematically explore these issues, an extension of the Standard Car Model widely used in research on ground vehicle platooning is proposed in this work. It represents the 3D double integrator kinematics of a UAS, further introducing linear first order input dynamics representing load factor tracking loops. Dynamic limitations of the UAS (stall angle, maximum thrust etc.) are represented by input magnitude saturations. This Standard Vehicle Model (SVM) is of low complexity compared to the more realistic model of a small UAS presented later, but is well suited to systematically study the impact of dynamic effects such as inner loop dynamics that pose serious challenges to the application of SMC techniques and in the development of new approaches to tackle these challenges. The equations of motion are given in a local North-East-Down tangent frame (index e). We denote at some points this frame as local inertial frame, since, as common in shortterm flight simulation, effects of the earths motion w.r.t. inertial space are neglected, as are local variations of gravity. The continuous-time vehicle position dynamics of vehicle i w.r.t. the local inertial frame follow from Newton's second law to

x e i (t) = (p e T i (t) v e T i (t)) T (1.1.1) ṗe i (t) = v e i (t) (1.1.2) ve i (t) = g e + R eg ,i n g i g (1.1.3)
where p e i (t) is the vehicle's position w.r.t. the NED frame, v e i (t) denotes its velocity, n g i are load factors due to aerodynamic forces including thrust acting on vehicle i expressed in the local guidance frame, R eg ,i performs the rotation from that frame to the NED frame, g e is the local gravity vector and g the norm of the nominal gravity vector converting load factors into accelerations.

Exogenous disturbances

The load factors acting on the aircraft are assumed to be composed of nominal load factor tracking outputs and unknown but bounded parasitic load factors due to atmospheric turbulence, wake vortices of neighboring UAS and imperfect tracking:

n g i = ñg i + n g w,i (1

Load factor tracking dynamics

It is assumed that inner load factor controllers track a decoupled first order dynamics reference model

ṅg i =    -τ -1 x 0 0 0 -τ -1 y 0 0 0 -τ -1 z    (ñ g i -sat(n g i (t), n g i (t), n g i (t))) (1.1.5)
where ng i are commanded load factors, ñg i is the state of the tracking loop, with the corresponding time constants τ p for p = x, y, z and n g i , n g i define time-varying box bounds on the admissible load factors. As will be shown later on, inner loop bandwidth and sampling time of the guidance laws pose the principal issues for continuous time sliding mode control. By varying these parameters, the SVM can cover both the simple double integrator kinematic model in its various forms used in the sliding mode literature (see e.g. [START_REF] Galzi | UAV formations control using high order sliding modes[END_REF]) and more realistic configurations. The specific vehicle dynamics are captured by the inner loop load factor controllers. Load factor tracking perturbations and tracking reference saturations n g i (t), n g i (t) are specific to a given vehicle and mission environment.

Electric glider UAS model

While the SVM is useful for preliminary studies, a more realistic 6DOF UAS model including rotatory dynamics, engine and actuator models is provided as a second major component of the FFB benchmark and is proposed to evaluate guidance and algorithms

Aerodynamics and actuators

There is an abundance of possibilities to model the aerodynamic coefficients of an aircraft, reaching from look-up tables to neural networks. Following the principle of minimum necessary complexity, here a simple model of the aerodynamic coefficients based on a truncated Taylor expansion and a quadratic drag polar has been adopted. The aerodynamic force and moment coefficients in the wind frame are given by

C D = C D 0 + (C L,0 + C L,α α) 2 1 πeΛ + C D δe δ e + C Dq q + C D δ f δ f (1.2.1) C Y = C y β β + C Yp p + C Yr r + C Y δa δ a + C Y δr δ r (1.2.2) C L = C L 0 + C Lα α + C Lq q + C L δe δ e + C L δ f δ f (1.2.3) C l = C l β β + C lp p + C lr r + C l δa δ a + C l δr δ r (1.2.4) C m = C m 0 + C mα α + C mq q + C m δe δ e + C m δ f δ f (1.2.5) C n = C n β β + C np p + C nr r + C n δa δ a + C n δr δ r (1.2.6)
Aerodynamic coefficients are adopted from [START_REF] Stolle | Towards vision-based autonomous cross-country soaring[END_REF], where the vortex lattice software AVL is used to obtain an aerodynamic model of a small electric glider airplane (Multiplex Cularis, depicted in flight in figure 1.2).

Control surface actuators and engine

A second order approximation of a typical miniature servo actuator has been determined in previous work [START_REF] Ameho | Du pilotage d'une famille de drones à celui d'un drone hybride via la commande adaptative[END_REF] by system identification:

δi δi = 0 1 -ω 2 n -2ζω n δ i δi + 0 ω n δ i,c (1.2.7)
with the natural frequency ω n = 62.8 s -1 and the damping coefficient ζ = 0.8. This model is augmented here by adding realistic saturations on deflection and deflection rate

δ i ∈ [- π 4 , π 4 ] rad (1.2.8) δi ∈ [-5.8, 5.8] rad s (1.2.9)
The engine is modeled as first order system

δen = τ -1 en (δ en,c -δ en ) (1.2.10)
where τ en = 0.27 -1 s -1 . The thrust T is then computed using a simple model based on momentum theory (from [START_REF] Beard | Small Unmanned Aircraft: Theory and Practice[END_REF]) as

T = ρ 2 S en C en ((k en δ en ) 2 -V 2 a ) (1.2.11)
with the propeller disc surface S en , the thrust constant C en and the engine constant k en .

Kinematics

Using 1.2.1-1.2.6, the aerodynamic forces and moments in the body frame are computed as

F a = R bw qS    C D C Y C L    +    T 0 0    (1.2.12) M a = R bw    qSb 0 0 0 qSc 0 0 0 qSb       C l C m C n    (1.2.13)
with the vector of forces F a and the vector of moments M a , assuming that the thrust vector is aligned with the body frame x axis and coincides with the center of gravity.

Translational accelerations in the body frame result to

vb = 1 m F a -(ω b × v b ) + R be g e (1.2.14)
where v is the translational velocity of the body frame w.r.t the inertial frame expressed in body frame axes, and ω b is the angular velocity of the body frame w.r.t the inertial frame expressed in body frame axes. The rotation rate of the earth can be neglected due to the short duration and range of the considered UAS flight maneuvers (the full benchmark takes about 4 minutes to complete at a nominal trajectory speed of 15 m/s) and our focus on local relative guidance. Angular accelerations about body frame axes follow to

ωb = -J -1 (ω b × (J ω b )) + J -1 M a (1.2.15)
with the inertia tensor J. With these accelerations, standard integration techniques and transformations can be used to propagate angular and translational velocities and positions over time to compute the aircraft's full state, which is given by

x = p e T v e T q T eb ω b T T (1.2.16)
with NED position p e , NED velocity v e , the attitude quaternion q eb and the vector of rotation rate about body axes ω b .

Wind disturbance modeling

The UAS is subject to exogenous disturbances due to atmospheric free air turbulence and the wake vortices of surrounding UAS. These are taken into account as additional system inputs, extending the input vector to

u = δ e δ a δ r δ en δ f v T w ω T w T (1.2.17)
with the control surface deflections δ e (elevator), δ a (differential aileron), δ r (rudder), δ en (throttle setting), δ f (flaps) and the induced wind v w and induced angular rates ω w .

Wind

The atmospheric wind is assumed to be a linear combination of a constant ambient part and a time-varying stochastic part, leading to induced velocities in the body frame

v b w (t) = R be (v e w,a + v e w,s (t)) (1.2.18)
and corresponding perturbations in airspeed V a and wind angles α, β and induced angular rates ω b w (t). Wind time series are generated according to the Dryden turbulence spectrum [START_REF] Moorhouse | US military specification MIL-F-8785C[END_REF].

Wake vortex disturbances A variety of approaches has been proposed to approximate the effects of trailing vortices on the following UAS, mostly based on modified Horseshoe Vortex models (HVM) (e.g. [START_REF] Hummel | Aerodynamic aspects of formation flight in birds[END_REF]), Vortex Lattice methods (VLM) (e.g. [START_REF] Saban | Simulation of Wake Vortex Effects for UAVs in Close Formation Flight[END_REF]), or even, very recently, real-time CFD approaches [START_REF] Tian | Wake Encounter Simulation and Flight Validation with UAV Close Formation Flight[END_REF]. For the proposed benchmark, a HVM with modified core model presented in ( [START_REF] Dogan | Modeling of aerodynamic coupling between aircraft in close proximity[END_REF]) is adopted. It is reported to provide predictions that are in good agreement with both VLM models and wind tunnel measurements while allowing for reasonable simulation times on a standard PC. From a set of induced velocities at discrete points on the airframe, by using an averaging scheme, an approximate induced airspeed vector and induced angular rates are generated in function of the inter-UAS separation vector, see [START_REF] Dogan | Modeling of aerodynamic coupling between aircraft in close proximity[END_REF] for details.

In figs. where P i , P j are the mean engine power in cruise flight of the follower and the predecessor, respectively. The lateral and vertical separations ∆y g o , ∆z g o are deviations from the position of minimum power draw.

Note that for the Cularis UAS model, the model predicts rather modest maximum gains of just over 6 %. This does not come unexpectedly due to the approximate nature of the vortex model and the quite complex aerodynamic interactions considered. Note, however, that the objective of the benchmark is not accurate quantitative prediction, but comparison of competing flight control strategies under identical conditions. To this end, a vortex interaction model that predicts the major trends is considered fully sufficient. Trim aileron deflections displayed in figures 1.3 and 1.4 suggest that in this particular setting, vortex disturbances are not a major factor for attitude control, since input resources are only reduced by maximally about 10 %.

Benchmark trajectory

To be feasible, a flight trajectory needs to be consistent with the dynamics of the UAS. Strictly, this implies continuity of flight path derivatives down to servo control inputs, taken even more strictly down to the second derivative of servo inputs, assuming second order servo dynamics. For a fixed wing UAS this translates into C 6 continuity. In reality, C 2 continuity, corresponding to continuous accelerations, is usually sufficient due to very small differences in resulting trajectories when requiring higher order continuity.

The proposed benchmark includes a trajectory composed of cubic Bézier splines, based on a two-step feasible path planning approach, presented by the authors of [START_REF] Yang | 3D smooth path planning for a UAV in cluttered natural environments[END_REF] for autonomous helicopter flight. After planning first a waypoint-based trajectory composed of straight lines, a series of C 2 continuous smoothing splines is computed to obtain a feasible approximation. Maximum curvature is a design parameter of the smoothing splines and allows to adapt it to specific vehicle constraints such as maximum lift.

The curvature of the benchmark trajectory is depicted in figure 1.5. The benchmark trajectory consists of a straight part to evaluate cruise flight performance, followed by subsequent climb and descent maneuvers, a 90°turn and, inspired by [START_REF] Yamasaki | Separate-channel integrated guidance and autopilot for a path-following UAV via high-order sliding modes[END_REF][START_REF] Looye | Helical flight path trajectories for autopilot evaluation[END_REF], an approximately helical part (see figure 1.6) to evaluate performance during simultaneous climb/descent and turn maneuvers. In the context of upwash exploitation, the cruise flight segment allows to compute useful figures of merit, since typically a major part of UAS long range missions is spent in cruise flight. The maneuvering segments of the trajectory add predecessor accelerations as additional disturbance and are suited for the study of maximum control errors as well as the propagation of maneuvering effects.

Arc length parameterization

Smoothing a waypoint-based trajectory following [START_REF] Yang | 3D smooth path planning for a UAV in cluttered natural environments[END_REF], the resulting curvature-limited sequence of cubic Bézier splines are parameterized with the spline parameter t ∈ [0, 1].

For evaluation of the formation position and velocity over time, an arc length parametrization is more convenient. Using the method proposed by the authors of [START_REF] Wang | Arc-length parameterized spline curves for real-time simulation[END_REF], we first approximate the arc length of each trajectory spline S i by integration of piecewise linear samples and then fit a piecewise polynomial cubic spline Si to discrete samples of the spline parameter over the arc length. For a given arc length the corresponding spline parameter is then found by evaluating Si .

Trajectory primitives

Based on the corner smoothing approach [START_REF] Yang | 3D smooth path planning for a UAV in cluttered natural environments[END_REF], a number of trajectory primitives have been implemented that can be chained to form arbitrary maneuver sequences in possible future extensions of the benchmark trajectory. Each trajectory is built up by adding primitives in an incremental manner, i.e. when adding a primitive, the heading of the last leg ψ and the end point P is modified by an certain increment, avoiding the need for keeping track of frame transformations. In the following we use the rotation matrix between the NED frame and a frame corresponding to an NED frame rotated about the NED z axis by a given heading angle defined by

R et (ψ) =    cos(ψ) -sin(ψ) 0 sin(ψ) cos(ψ) 0 0 0 1s    (1.
∆P = R et ( ψ + ∆ψ)    d 0 0    (1.3.2)
Climb-Descent The climb primitive is based on adding a waypoint to the trajectory in function of a desired climb angle γ, altitude increment δh and trajectory heading increment ∆ψ. The waypoint increment results then to

∆P = R et ( ψ + ∆ψ)    ∆h tan γ 0 -∆h    (1.3.3)
Turn A turn adds two successive waypoint increments according to a desired turn angle ∆ψ and turn leg distance d as

∆P 1 = R et ( ψ)    d 0 0    , ∆P 2 = R et ( ψ + ∆ψ)    d 0 0    (1.3.4)

Helix

The helix maneuver primitive is parameterized by turn radius r, climb or descent altitude per turn ∆h, vertical center unit vector c indicating sense of rotation and number of turns n. Due to their non-constant curvature, Bézier splines can only approximate a circular helix. The approximation can be made arbitrarily close to the true helix by increasing the number of waypoints per turn N , with the corresponding increase in memory footprint. Starting from the initial heading ψ the initial radius vector and the initial waypoint increment are computed to

r 0 = -rc × R et ( ψ)    1 0 0    (1.3.5) ∆P 1 = R et ( π N )r 0 -r 0 (1.3.6)
and the sequence of waypoint increments is given recursively by

∆P i = R et ( 2π N )∆P i-1 +    0 0 -∆h N    for i = 2...nN (1.3.7)

Visualization

A lightweight visualization environment based on the Java jmonkey game engine [START_REF]jmonkeyengine.org Java Game Engine[END_REF] has been developed to display the attitude and positions of multiple UAS in a synthetic 3D environment. It has proven to enhance productivity while developing and debugging control laws as well as the UAS simulation model itself. It is interfaced with the Simulink dynamics simulation via a UDP (User Datagram Protocol) network link. Since the visuals have to be updated only at a moderate rate of about 30 Hz for human perception, sending UDP packets only adds a small overhead to the overall simulation time.

The tool currently provides several UAS 3D models (glider, electric glider, quadcopter) as well as the possibility to display wake vortex cores.

Computational effort Formations of three UAS have been simulated on an Intel i7 4-core machine close to real time in Simulink ® normal simulation mode. Running the benchmark in Simulink ® Accelerator mode, where a standalone executable is compiled for the host machine, allows for considerable speedups. It comes however with the downside that control laws have to be implemented in the limited subset of Simulink ® /Matlab ® compatible with code generation.

On guidance and control algorithms for tight formation flight

A variety of guidance architectures has been considered in the formation flight and vehicles platooning community, from fully centralized ones over those employing leader information to ensure string stability, circular architectures and bidirectional to fully decentralized ones relying only on relative states between neighboring vehicles. In this work we consider a unidirectional predecessor-follower configuration, as it minimizes data exchange within the formation and is convenient for visual localization.

Main objectives

Tight formation flight confronts us with two principal control objectives: avoiding collisions entirely and maximizing the aerodynamic gain from the neighbor's wake. For the first objective, obtaining a useful, i.e. not overly conservative upper bound on the maximum control error is desirable.

For the second objective it is convenient to consider the loss in upwash gains introduced by imperfect tracking of the sweet spot or equivalently the added power consumption. Loss in energy savings correspond to integrated loss in propulsion power savings, which are proportional to relative position tracking errors. The relation can locally be approximated as quadratic (see figs. 1.3,1.4)

∆P k ≈ (c y ∆y g 2 o,k + ∆z g 2 o,k ) (2.1.1)
The added energy consumption for a configuration λ over a flight of N samples with sample time T s then being approximated by

∆E λ ≈ T s N k=1 ∆P λ,k (2.1.2)
we can write the relative efficiency of two configurations λ and κ as

∆e ≈ T s N k=1 ∆P λ,k T s N k=1 ∆P κ,k (2.1.3) = 1 N N k=1 (c y ∆y g 2 o,λ,k + c z ∆z g 2 o,λ,k ) 1 N N k=1 (c y ∆y g 2 o,κ,k + c z ∆z g 2 o,κ,k ) (2.1.4)
The weighted mean squared (WMS) tracking error

e λ (c y , c z ) = 1 N N k=1 (c y ∆y g 2 o,k + c z ∆z g 2 o,k ) (2.1.5)
is therefore a simple and meaningful performance index to compare the energy efficiency of different configurations (such as different guidance laws). Recall that in the longitudinal channel, larger errors can be tolerated since the upwash gain is much less sensitive than it is w.r.t. lateral and vertical variations. The predecessor tracking problem we deal with in this chapter is thus the following:

Minimize the weighted sum of squared tracking errors given by relation (2.1.5), while providing guaranteed, i.e. deterministic bounds on the maximum deviation from the nominal relative position.

This tracking problem is generally harder for smaller UAS than larger aircraft as we elaborate on in the next section in some more detail.

Scaling effects affecting small UAS

Tight formation flight demonstrating large upwash gains under automatic control has been demonstrated on a pair of supersonic F/A-18 fighter aircraft in 2001 [START_REF] Jake | F/A-18 Performance Benefits Measured During the Autonomous Formation Flight Project[END_REF]. Rare demonstrations with small UAS that are close to TFF tracking requirements have been reported more than ten years later [START_REF] Wilson | Vision-aided guidance and navigation for close formation flight[END_REF], raising the question of fundamental obstacles, in more informal terms: if it has been done on large aircraft, why not on small ones? Scale analysis allows some order-of-magnitude insight into how airframe size impacts the TFF problem.

Scaling of sensitivity towards turbulence

Atmospheric turbulence and wake vortices are the predominant source of disturbance in tight formation flight. Disturbances in the vertical translational dynamics due to atmospheric turbulence can be separated into

• disturbances of direction of the airspeed vector, modifying the angle of attack with the corresponding modification of aerodynamic forces and moments

• disturbances of the airspeed norm

Considering the airspeed vector expressed in the aircraft body frame, the angle of attack is given by

α = sin -1 w V a (2.1.6)
where V a = (u v w) T . Since α is limited to small angles α < α stall , equation (2.1.6) can be approximated by

α ≈ w V a (2.1.7)
Adopting a truncated first order Taylor expansion of the lift coefficient, the vertical load factor due to lift, using 2.1.7 is approximated by

n z ≈ (C L 0 + C Lα α)Sb ρ 2 V 2 a gm (2.1.8) ≈ (C L 0 V a + C Lα w)Sb ρ 2 V a gm (2.1.9)
The partial derivative of 2.1.9 w.r.

δn z δV a = 1 V a (2.1.13)
It is therefore a useful and simple comparative metric. Typical flight conditions of the small UAS considered in the FFB benchmark (see chapter 1) and the NASA AFF program [START_REF] Jake | F/A-18 Performance Benefits Measured During the Autonomous Formation Flight Project[END_REF][START_REF] Bowers | An Overview of the NASA F-18 High Alpha Research Vehicle[END_REF] and the corresponding sensitivities are given in tab. 2.1. Note that for the Cularis UAS, an airspeed perturbation of the same magnitude causes load factor perturbations that are roughly ten times as high as those acting on the F/A-18 aircraft. Similar arguments can be made for the lateral channel and the rotatory dynamics. This snapshot result sheds some light on the relative difficulty of performing TFF with small UAS.

A brief state of the art

Despite a considerable interest in tight formation flight strategies for small UAS, only a few flight test results are reported so far. This is in contrast to a larger body of published experimental results on automatic TFF with manned aircraft. After early advancements [START_REF] Park | A new nonlinear guidance logic for trajectory tracking[END_REF] achieving meter-level formation flight of a heterogeneous pair of two small UAS, incremental improvement could be gained in the following years, see e.g. [START_REF] Gu | Design and flight testing evaluation of formation control laws[END_REF][START_REF] Rice | Autonomous close formation flight control with fixed wing and quadrotor test beds[END_REF], and a broad variety of guidance schemes tailored to predecessor following has been developed.

In an attempt at automatic aerial drogue docking with two identical small UAS (b = 2 m), the authors of [START_REF] Wilson | Vision-aided guidance and navigation for close formation flight[END_REF] report impressive achieved position errors of (∆x T RM S ∆y T RM S ∆z T RM S ) T = (0.42b 0.44b 0.22b) T in cruise flight under an average wind of 30 % of the airspeed, using proportional guidance laws augmented with a feed forward scheme, operating on inertial speed as interface variable. The use of monocular machine vision in combination with single-frequency GNSS makes the presented ap-proach well suited for low-cost UAS. Both [START_REF] Park | A new nonlinear guidance logic for trajectory tracking[END_REF] and [START_REF] Wilson | Vision-aided guidance and navigation for close formation flight[END_REF] choose velocity as interface variable between inner and outer loops In all cases collision avoidance guarantees are not taken into account explicitly and control synthesis involves manual gain tuning. These experimental results suggest that performance close to the requirements of meaningful upwash exploitation and aerial docking is available based on existing guidance and control techniques is within grasp. However, improvement is still necessary, achievable for instance by better system identification and careful control gain tuning.

On the other hand, sliding mode control in theory promises excellent performance and disturbance robustness while requiring little knowledge about the controlled system. Sliding mode control gains can systematically be derived from disturbance bounds. As is well known, adapting sliding mode to the constraints posed by real-world system invariably deteriorates the excellent performance of ideal sliding mode. We consider it worthwhile, as a research effort complementary to established guidance techniques to explore techniques to minimize the extent of this performance deterioration. Excellent position tracking performance could be the reward, as our simulations suggest, see section 2.6. In the following, we first propose a baseline control architecture -adapted to the hardware resources of small UAS -accompanied by baseline LQ guidance and control laws. We then proceed to the design of a novel predictive discrete sliding mode control law in section 2.5.1 that can considerably reduce position tracking errors -find simulation results in section 2.6 -and requires less tuning effort from the control designer than the LQ guidance laws.

Notation and terminology

In slight contrast to standard terminology, we denote the separation between UAS along the x axis of the guidance frame as longitudinal, along the y axis as lateral, and along the z axis as vertical separation. This is due to the aforementioned less critical importance of downstream inter-vehicle separation.

Guidance and control architecture

Maintaining a position in the wake of a preceding aircraft close to that of maximum upwash energy gain constitutes a 4D trajectory tracking problem. Not only a position trajectory needs to be followed (the path tracking problem), but each position needs to be visited at a specific time. By posing the problem in a moving frame attached to the UAS's nominal position, or, equivalently, considering the deviation from a possibly time varying reference position, the problem can, without loss of generality, be reduced to a 3D position regulation problem subject to a number of disturbances the control laws need to reject. Tracking trajectories within the formation, e.g. during a leader switch in a rotating Figure 2.1: Guidance and control architecture of follower i and a predecessor j leader configuration, can be handled in the same way by treating the change in reference position as further disturbance. We make no assumptions about the orientation of the formation trajectory w.r.t. to gravity other than it is feasible. We make some assumptions on what constitutes the available control inputs of a typical small UAS.

Control resources: assumptions

The minimally available control inputs of a typical fixed-wing UAS configurations are thrust setting and commanded elevator, aileron and sometimes rudder deflection. These inputs naturally limit the control bandwidth (we use the term informally here, since it only applies to linear systems in its strict sense) by which aerodynamic translational forces can be generated and disturbances can be counteracted, since lift is influenced by rotating the airframe about the y axis, and generating side forces involve either banking or yawing. The closed loop rotatory dynamics obviously depend on airframe parameters and available servo actuators, but in very general terms, both longitudinally 1 , laterally and vertically control bandwidths of the same order of magnitude are available. Considerable bandwidth in the critical lateral and vertical channel can be gained by adding control surfaces that have a direct impact on aerodynamic lift and side force, rendering the aircraft fully actuated. This approach has been taken by [START_REF] Park | Avionics and Control System Development for Mid-Air Rendezvous of Two Unmanned Aerial Vehicles[END_REF]. Side force generators add the benefit of further decoupling lateral and vertical control. There is thus a trade off between available control bandwidth and the design effort of retrofitting direct force generating surfaces. This effort can be considered significant especially in the lateral channel, adding furthermore additional drag, whereas flaps for Direct Lift Control (DLC) can be added comparatively simple. We therefore assume in this work that flap control surfaces are available. Note that some flying wing configurations of low wing sweep do not permit using flaps for DLC.

Two-layer nested architecture

Although by now means imposed by aircraft dynamics it is widely common in flight control to adopt a multiple-loop architecture, subdividing the problem into "guidance" and "control" problems. Arguments in favor of nested control loops are of practical nature. Both control design and implementation are facilitated, since inner loops can be designed and tested separately. Inner loop control laws can be re-used and outer guidance loops can be developed on simplified models of the inner closed-loop dynamics. Nested architectures provide an abstraction layer bearing strong resemblance to abstraction layers found in software design and share the same benefits, i.e. reducing the complexity of a system to manageable proportions. As a further benefit of multiple loops, guidance laws can be implemented at the typically moderate observation rates of localization sensors, reducing the need for high-rate state estimation. A great many existing UAS trajectory or path tracking control architectures are based on nested control laws, often relying on timescale separation arguments. As is argued in this chapter, by the right choice of separation into outer and inner loop states, full use can be made of the aforementioned advantages by reducing the inner loop closed loop dynamics to an uncertain, but purely kinematic system. A popular choice of interface variable2 are Euler angles, since in planar cruise flight the pitch angle Θ is approximately related to the angle of attack α and thus aerodynamic lift; the bank angle is approximately linearly proportional to lateral load factors (see section 2.4.3) for small angles. During trajectory transients, e.g. climb maneuvers, this Θ ≈ α assumption holds less, invariably deteriorating control performance. Tracking angle of attack at the inner loop level circumvents this issue, but requires additional instrumentation. Furthermore, from the perspective of the outer loops, the inner loop appears as time-varying system depending on a number of aerodynamic coefficients such as C Lα of which usually only approximations from system identification are known. Since both inner and outer loops are affected by uncertain aerodynamic coefficients, robust control techniques that can deal with these parametric uncertainties need to be employed on both levels. It would appear desirable to choose a set of interface variables that contain model uncertainty at the inner loop level. Both load factor and velocity are candidates for such a kinematic inner/outer loop interface. Estimating inertial velocity onboard small UAS relies however exclusively on GNSS Doppler and carrier phase observations. A temporary loss of GNSS observations renders the inner loops nonoperational. Load factors, on the other hand, can be estimated under loss of GNSS, although with reduced accuracy due to increasing AHRS errors without GNSS. 3D load factors as interface variables, although employed e.g. by [START_REF] Peddle | Acceleration based manoeuvre flight control system for unmanned aerial vehicles[END_REF], appear to have found negligible attention in the flight control community, while providing consider-able advantages that are exploited in this work. Nonlinear dynamic inversion (NDI) based inner loops such as proposed in [START_REF] Gu | Design and flight testing evaluation of formation control laws[END_REF] are a similar approach, making however the strong assumption that aerodynamic forces can be generated instantaneously, i.e. that the vehicle is fully actuated. Load factor tracking, on the other hand does not preclude taking tracking dynamics -due to rotatory and actuator dynamics -explicitly into account. Being able to track load factors allows to base the design of outer laws on a purely kinematic design model. In the ideal case, the dynamics of the vehicle are completely governed by the inner loops. Note that uncertainties of the aerodynamic model of the vehicle affect the outer loops only insofar as they might lead to inner loop tracking errors. A load factor interface is vehicle agnostic, insofar as it can accommodate fixed wing as well as helicopter and multirotor configurations and could even be used for underwater vehicles. The specificities of the different vehicles manifest themselves in the achievable closed loop load factor tracking dynamics and time-varying constraints on admissible load factors. Guidance laws can be developed on simple kinematic models, provided accurate load factor tracking can be implemented on the concerned vehicle. Furthermore an existing interface based on Euler angles and thrust can be cast as a load factor interface by simple transformations. Note that the above arguments do not constitute a claim that load factor tracking is an inherently superior choice. We claim however that it is a convenient and elegant choice, largely reiterating arguments given in [START_REF] Peddle | Acceleration based manoeuvre flight control system for unmanned aerial vehicles[END_REF], since it allows to confine aircraft-specific dynamics to the inner loops, thus providing a very useful abstraction layer for guidance loop design. In the following we consider a two layer guidance and control architecture as depicted in figure 2.1. Its main elements are load factor tracking loops exposing commanded load factors as virtual inputs and a novel predictive discrete sliding mode guidance law complemented by a disturbance estimation scheme. Load factor control as well as guidance laws operate in a local guidance frame defined in the sequel.

Baseline guidance and control laws 2.4.1 Guidance frame

The system to be controlled is a formation of n UAS flying in an arbitrary pattern. A variety of guidance frames has been proposed in the literature, ranging from a planar frame aligned with the follower's velocity ( [START_REF] Gu | Design and flight testing evaluation of formation control laws[END_REF]) to the predecessor's body frame ( [START_REF] Schumacher | Adaptive control of uavs in close-coupled formation flight[END_REF]). The primary objective of tight formation flight for energy saving purposes is to keep each UAS at the position of maximum efficiency in the wake of its predecessor. Being induced by the aerodynamic flow, the wake vortices are approximately aligned with the predecessor's instantaneous wind frame, plus the effect of trajectory curvature. For

x b,j z b,j x g,j z g,j z b,i ∆p c x e z e ∆p x b,i

Nominal predecessor vortex frame

Since small UAS typically are not equipped with sensors for angle of attack and side slip angle, a local guidance frame (index g) is introduced that approximates the predecessor's vortex frame. Its x axis unit vector u x,g is aligned with each predecessor's NED velocity vector. Its y axis should be aligned with the y axis of the predecessor's body frame, since the central axis of the vortices lies somewhere on the predecessor's wing. To avoid having to communicate the NED speed and attitude of each UAS to its follower, the nominal speed and acceleration vector of each UAS on the current formation trajectory is used to derive a nominal guidance frame. These frames are computed sequentially for each UAS, since the attitude of each nominal local guidance frame defines the nominal NED position of the next one. For two subsequent UAS, predecessor p and follower f , the nominal NED speed of f is computed as

v e f = v e p + ω e p × ∆p e f (2.4.1)
with the commanded separation vector ∆p e f . Acceleration is given by

ve f = ve p + ω e p × (ω e p × ∆p e ) (2.4.2)
and the instantaneous rotation rate ω e f of the guidance frame due to trajectory curvature by

|ω e f | = | ve f | |∆p e | 1 2
(2.4. 3)

ω e f = |ω e f | ve f × v e f | ve f × v e f | (2.4.4)
The x axis unit vector of the local guidance frame results then to

u x = |v e f | -1 v e f (2.4.5)
The orientation of the z axis unit vector u z is derived from the simplifying assumption that gravitational acceleration and centrifugal forces due to trajectory curvature are compensated for by the thrust and drag along x g and by the aerodynamic lift in a plane normal to x g . Further assuming that the aerodynamic lift Z lies in the x-z plane of the body frames, first the total centrifugal and gravitational acceleration acting on the aircraft is computed as

a t = |ω p | 2 r + (0 0 g) T (2.4.6)
The z unit vector u z is then computed by projecting and normalizing the total acceleration on the z-y plane of the g frame

a t,zy = a t - a t • u x u x • u x • u x (2.4.7) u z = |a t,zy | -1 a t,zy (2.4.8)
and u y completes a right-handed Cartesian coordinate frame

u y = u z × u x (2.4.9)
leading to the corresponding rotation matrix

R eg = u x u y u z (2.4.10)
composed of the unit axes of the guidance frame, see figure 2.2.

Differential dynamics of a pair of vehicles

It complicates modeling, but simplifies control allocation to write vehicle dynamics in the guidance frame. The vehicle position in the guidance frame is given by

p g i (t) = R ge (t)(p e i -p e eg ) (2.4.11)
where p e eg is the position of the g frame w.r.t. the e frame expressed in the e frame. The relative position error of vehicle i + 1 w.r.t. its predecessor i is then given by

∆p g i (t) = R ge (t)(p e i+1 -p e eg ) -R ge (t)(p e i -p e eg ) -∆p g c,i (t) (2.4.12) = R ge (t)(p e i+1 -p e i ) -∆p g c,i (t) (2.4.13)
with its first derivative

∆ ṗg i (t) = R ge (t)(v e i+1 -v e i ) -Ω g (t)∆p g -∆ ṗg c,i (t) (2.4.14)
and second derivative

∆p g i (t) =R ge (t)( ve i+1 -ve i ) + Ṙge (t)(v e i+1 -v e i ) -Ω g (t)∆ ṗg -Ωg (t)∆p g -∆p g c,i (t) (2.4.15)
where ∆p g c and its first two derivatives define the desired trajectory of the follower in the guidance frame and where Ω g is the cross product matrix ( [START_REF] Stevens | Aircraft control and simulation: dynamics, controls design, and autonomous systems[END_REF]) of the angular velocity ω g = ω g ,1 ω g ,2 ω g ,3

T of the guidance frame w.r.t. the inertial frame

Ω g =    0 -ω g ,3 ω g ,2 ω g ,3 0 -ω g ,1 -ω g ,2 ω g ,1 0    (2.4.16)
and using

Ṙge = -Ω g R ge (2.4.17)
Introducing

∆v g i (t) =R ge (t)(v e i+1 -v e i ) -∆ ṗg c,i (t) (2.4.18) ξ v =Ω g (t)∆p g (2.4.19) ξ a =R ge (t) ve i -Ω g (t)(∆v g i + ∆ ṗg c,i (t)) -Ω g (t)∆ ṗg -Ωg (t)∆p g -∆p g c,i (t) + g g + gn g w,i+1
(2.4.20)

and rearranging, the assembled dynamics in the g frame can be written in a linear form convenient for control design as

∆ ṗg i (t) = ∆v g i (t) + ξ v ∆ vg i (t) = gn g i+1 + ξ a (2.4.21)
with the control input n g i+1 and the partly known disturbances ξ v and ξ a . The explicit dependence on time is dropped in the following for better readability. Four different sources of disturbances become apparent here:

• disturbances corresponding to inevitable load factor tracking errors of the follower w.r.t. inertial space (n g w,i+1 )

• disturbances corresponding to the rotating guidance frame (i.e. Ωg (t)∆p g )

• disturbances corresponding to predecessor accelerations ( ve i )

• disturbances corresponding to modifications of the commanded relative position (∆p g c,i (t))

It is the objective of the guidance laws to track desired positions w.r.t. the predecessor UAS. In the following we first present a set of baseline control laws based on the Linear Quadratic Regulator (LQR) for the system given by equation (2.4.21) and then turn our attention to potentially more powerful sliding mode guidance techniques.

LQ baseline guidance laws

To provide a standard to which to compare new control approaches, we design a set of baseline LQ control laws as part of the FFB benchmark. For controller synthesis, the respective subsystems are extracted from the full linear model. The structure of the control laws is based on the usual assumptions of weak coupling between lateral and longitudinal dynamics and timescale separation of attitude dynamics and translational dynamics.

The synthesis procedure is the same for all four subsystems. It is shortly recalled here, and the specific design systems for every controller are given in the respective subsections. The time-invariant continuous time linear quadratic regulator minimizes the following quadratic performance index

J = ∞ t0 (x T Qx + u T Ru)dt (2.4.22)
with the positive semi-definite state cost matrix Q and the positive definite input cost matrix R for a given linear system

ẋ(t) = Ax(t) + Bu(t) (2.4.23)
To track a given reference state x c , the tracking error dynamics are derived (dropping in the following the explicit dependence on time) as

ẋe = ẋ -ẋc (2.4.24) = Ax + Bu -ẋc (2.4.25)
The reference state is either constant, as in the static relative position tracking case, or generated by the outer loop controller (e.g. the commanded bank angle). Its derivative is here simply treated as unknown disturbance and set to zero for controller design, leading to the state error dynamics

ẋe = Ax + Bu (2.4.26)
that are identical with the state dynamics. Tracking is thus achieved by a simple change of state coordinates

x e = x -x c (2.4.27)
To cope with steady state tracking errors, integral action is added by augmenting the system (2.4.27) with integral states

ẋe ẋi = A 0 G 0 x e x i + B 0 u (2.4.28)
where G is a sparse matrix defining the states that are selected for integral action. As is common practice, wind disturbances and actuator and engine dynamics are neglected for the LQR design systems. A linear model including models of actuators and engine is used for fast evaluation of candidate controllers during controller tuning. Note that states and inputs are deviations from the operating point where the design model has been obtained by linearization and that for implementation the corresponding operating point trim values need to be added.

Predecessor tracking

Longitudinal The control inputs of the longitudinal position control loop are the commanded longitudinal load factors n x and n z . In the absence of external disturbances the integral LQ baseline controllers regulate the relative position and velocity errors asymptotically to zero. The design system in the guidance frame is given by

          ∆ ṗg x ∆ ṗg z ∆ vg x ∆ vg z ∆ ṗg x,int ∆ ṗg z,int           =    0 2 I 2 0 2 0 2 0 2 0 2 I 2 0 2 0 2              ∆p g x ∆p g z ∆v g x ∆v g z ∆p g x,int ∆p g z,int           +    0 2 gI 2 0 2    n x n z (2.4.29)
where I p denotes a unity matrix of size p × p and 0 p denotes a matrix of zeros of size p × p. Note that virtual accelerations due to the rotating guidance frame are lumped into the unknown disturbances here and neglected in the design model.

Lateral

The lateral control strategy is derived from the coordinated turn assumption

F e a,z ≈ mg (2.4.30
)

F e a ≈ mg cos Φ (2.4.31) ⇒ n g y,c ≈ sin Φ g cos Φ (2.4.32) ⇒ Φ g c = sin -1 (n g y cos Φ) (2.4.33)
thus lateral load factors in the local guidance frame are generated by inclining the approximate lift vector F e a via the bank angle w.r.t. the guidance frame, Φ g . The resulting simple LQ design system is given by

   ∆ ṗy ∆ vy ∆ ṗy,int    =    0 1 0 0 0 0 1 0 1       ∆p y ∆v y ∆p y,int    +    0 g 0    n y,c (2.4.34)
and the commanded bank angle in the local guidance frame is computed using equation (2.4.33).

We proceed to design a set of integral LQ load factor controllers for the vertical and lateral channel.

Some elements of load factor tracking loop design

The first set of controllers tracks pitch and bank angle tracking commands, representing a rather classic approach, and then applies transformations to approximate vertical and lateral load factor tracking. It is part of the benchmark presented in chapter 1. In a second step, an integral LQ control for vertical load factors is presented based on an augmented state and employing Direct Lift Control for improved tracking performance.

LQ baseline load factor tracking

Commercial as well as open-source UAS autopilots typically provide some kind of basic attitude tracking laws e.g. based on PID control. We keep this assumption for the baseline control laws, emulating the presence of standard inner loop control laws and design a set of LQ attitude tracking laws, but apply transformations to expose a load factor tracking interface to the guidance laws.

Attitude tracking

The inner loops consist of separate pitch angle and bank angle tracking laws. Both angles are not the usual Euler angles w.r.t. to the NED frame but are the instantaneous rotations about the y and x axes of the body frame required to rotate the body frame into the desired attitude. Assuming that the derivatives of the commanded attitude are sufficiently bounded, the body frame stays close to its desired attitude and issues arising at Euler angle singularities are avoided. That being said, the commanded attitude rotation matrix w.r.t. to the NED frame is composed as

R eb,c = R eg R bg T ,c (Φ g c , Θ g c ) (2.4.35)
where

R bg T ,c (Φ c , Θ g c ) =    cos Φ g c sin Φ g c cos Θ g c sin Φ g c cos Θ g c -sin Φ g c cos Φ g c cos Θ g c cos Θ g c sin Θ g c 0 -sin Θ g c cos Φ g c    (2.4.36)
Note that the heading angle in the guidance frame is not actively tracked and is set to zero. Virtual angular velocities due to the rotating guidance frame are lumped into the unknown disturbance and neglected in the design models.

Longitudinal The pitch attitude design system is given by

   Θg q Θg int    = A lon,att 0 1 0 0    Θ g q Θ g int    + B lon,att δ e (2.

4.37)

Vertical load factor tracking Approximate nominal vertical load factor tracking can be achieved by assuming Θ g ≈ α in level cruise flight and using the nominal lift coeffi-cient 1.2.3 to compute Θ g c , dropping small coefficients to

Θ g c = C -1 Lα ( n z,c g qSb -C L 0 ) (2.4.38)
Longitudinal load factor tracking Approximate load factor tracking in x is obtained by inversion of the engine model 1.2.11, i.e. solving equation (1.2.11) for the engine command δ en,c given a commanded thrust force T c = n x,c gm.

Lateral load factor tracking

The lateral attitude law fulfills two requirements: tracking a commanded local bank angle and adding damping to the weakly damped dutch roll mode. The second objective is achieved by adding the yaw rate r to the state vector of the design system and choosing an appropriate weight in the LQR synthesis matrix Q. The design system is thus given by:

     ∆ Φ ṗ ṙ ∆ Φint      = A lat,att 0
1 0 0 0      ∆Φ p r ∆Φ int      + B lat,att δ a δ r (2.4.39)

Incorporating Direct Lift Control

Direct Lift control (DLC) has been considered as early as the sixties [START_REF] Klein | A preliminary investigation of the use of lift control in maneuvering flight[END_REF] to improve fixed wing aircraft responsiveness in the vertical channel. Used in manned refueling, the authors of [START_REF] Mcneill | A flight study of the use of direct-lift-control flaps to improve station keeping during in-flight refueling[END_REF] report improvements in station holding by flaperons. In [START_REF] Lombaerts | Design and flight testing of nonlinear autoflight control laws incorporating direct lift control[END_REF] Direct Lift Control by dedicated control surfaces is used for improved flight path following.

In another very recent contribution [START_REF] Jaiswal | Adaptive Longitudinal Control of UAVs with Direct Lift Control[END_REF] an adaptive NDI scheme is successfully employed using DLC to improve altitude tracking of a UAS under the assumption of separate flaps. Using DLC, aerodynamic lift can be influenced with considerably higher bandwidth than through the pitch channel. The Cularis UAS model that is part of the FFB benchmark features separate flaps for improved glide performance while thermaling and shortened landing distances. Here we exploit these control surfaces for Direct Lift Control. The most simple control structure consists of tracking vertical load factor commands solely by DLC, while stabilizing the short period motion by some auxiliary control. DLC lift resources are however limited due to mechanically limited flap deflections and flap size. The two input channels, pitching -slower, but with maximum lift only limited by the stall angle and available thrust -and DLC surfaces -fast but of limited magnitude -have thus complementary properties when it comes to lift control. As we will see in this section, by augmenting the linearized load factor dynamics with auxiliary states, LQ state feedback control can be designed that exploits these complementary properties in a systematic way via the usual LQR tuning matrices. Consider the short period dynamics captured by the vertical load factor n z and the pitch rate q [50]

ṅz q = A n z q + B δ e δ f (2.4.40) n z q = C n z q + D δ e δ f (2.4.41)
where we assume that C = I 2 , i.e. direct observations of load factor and pitch rate are available from the onboard AHRS. Note that the feedthrough matrix D has nonzero entries, since both control entries generate direct lift. We extend the state vector by defining the input derivatives as new virtual control inputs, rendering D zero. We further add the integral of n z to the states for zero steady steady state tracking error. Furthermore, we would like the control to employ δ f only for the high-frequency portion of the tracking task, to avoid saturation. This is achieved by including the integral of the flap deflection δ f as additional state, leading to the augmented LQR design model

          ṅz q δe δf n z δ f           =        A B 0 2 0 2 0 2 0 2 1 0 0 1×2 0 1×2 0 1×2 0 1 0 1×2                  n z q δ e δ f n z δ f           +    0 2 I 2 0 2    δde δf (2.4.42)
The LQR cost-to-go matrices Q, R are selected as diagonal matrices. The diagonal entries of the tuning matrices have an intuitive appeal, as most have a direct impact of a specific characteristic of the resulting control, such as an increase in Q 2,2 increases the short period damping etc. See figure 2.3 for a tracking example. Note that due to the integral penalization of flap deflections δ f , flaps are only employed at the beginning of a tracking move and rapidly return to zero, while the lower frequency portion of the load factor is carried by aerodynamic lift due to angle of attack, indicated by non-zero steady state δ e .

Benchmark performance

When applied to the FFB benchmark, these basic guidance and control laws provide good performance under cruise flight conditions, but insufficient disturbance rejection during predecessor maneuvering, see simulation results in section 2. [START_REF] Beukenberg | Aerodynamics, performance and control of airplanes in formation flight[END_REF].

In what follows, we review the applicability of sliding mode control to the relative guidance problem. We present predictive discrete sliding mode guidance laws that provide improved tracking performance both in cruise flight and during maneuvers.

Predictive discrete sliding mode guidance laws

In most general terms, it is the objective of control to enforce desirable system dynamics. Sliding Mode Control (SMC) approaches this problem by first designing a hypersurface -the sliding surface (SSF) -in state space that represents desirable, stable closed-loop dynamics of the type

0 = σ(x, x (1) , ...x (n) ) (2.5.1)
where x (p) = d p x dt p and n is the highest derivative of x involved in the definition of the desired dynamics. This general definition is by now means exhaustive and often output based sliding surfaces are employed. Here we assume that estimates of the relative position and velocity states are available and focus on state-based sliding surfaces. To give a more concrete example, consider a perturbed first order system ẋ = u + w and a linear sliding surface

0 = λx + ẋ (2.5.2)
The general disturbances w can lump together nonlinearities, the effects of modeling errors as well as exogenous disturbances. This prescribed behavior forms a hypersurface (in this simple example a line) in state space. In a control synthesis step, a control in-put u is designed that is guaranteed to drive the system states towards this surface even under worst-case disturbances w (reaching phase), mostly using Lyapunov approaches. Once reached, the control keeps the system states confined to the SSF (sliding phase). While moving along the sliding surface, the closed loop system is indistinguishable from the design system represented by the sliding surface. In our example, the system evolves according to

ẋ = -λx (2.5.3)
which is stable for λ > 0. While greatly differing in detail, all sliding mode approaches are based on these basic principles. Note that the underlying idea of SMC is quite different from other control techniques and a quite audacious one: to instantly eliminate disturbances acting on a system by a counteracting control input. In ideal sliding mode, guaranteeing system motion towards the sliding surface at all times implies infinitely fast switching on the sliding surface of a control input of sufficient magnitude to overpower the worst possible disturbance driving the system away from the sliding surface.

It is this worst-case paradigm in concert with the use of continuous-time design models that causes at once the great robustness properties of sliding mode control for certain select classes of dynamic systems with matched3 uncertainty and potentially very poor performance for many classes of real-world aerospace systems. For an illustrative example, consider control of the vertical position of an aircraft. Disturbances enter as accelerations. In order to generate counteracting accelerations, control inputs have to rotate the aircraft about the y axis, adding at least two integrators between disturbances and inputs. Even in the presence of DLC devices, actuator dynamics render disturbances unmatched. What is more, guidance laws are implemented at finite discrete sampling rates, invalidating the assumption of instantaneous input switching. Both effects, as is well known, lead to chattering, high frequency oscillations about the sliding surface and in the control inputs.

The approximate nature of both its basic assumptions -matched disturbances and instantaneous switching -on most physical mechanical systems poses the major obstacle to implementation of continuous time sliding mode control and a number of strategies have emerged to improve on this fundamental issue. Boundary layer smoothing approaches are based on replacing the discontinuous control inside a boundary layer around the sliding surface by some linear or nonlinear continuous control law. Continuous approximations of the signum function such as the saturation or the sigmoid function are popular examples of this class of chattering attenuation techniques. While often successfully applied in practice, an important drawback results from the fact that ideal sliding mode is lost even in the disturbance free continuous time case.

A more recent approach is based on including higher-order derivatives of the sliding variable into the feedback. Higher order sliding mode control of order n forces not only the sliding variable, but its time derivatives up to the order of n -1 to zero. Input discontinuities appear in the (n -1)th derivative of the control input. The necessity for higher-order derivatives of the sliding variables for feedback purposes, introducing the difficult task of real time signal differentiation, poses the principal implementation issue of HSOM techniques. HOSM observers have been proposed for this purpose, which however remains a challenging task in the presence of discrete sampling and noise. As a notable exception for systems of relative degree one, the Super-Twisting controller (STSMC) does not require derivatives of the sliding variable to enforce second order sliding mode, i.e.

σ = σ = 0 (2.5.4)
in finite time, making it a very popular choice in SMC applications. It is often proposed as a remedy to recover perfect second order tracking in the ideal case, while providing continuous implementable control input signals.

Discrete Sliding Mode Control (DSMC), treated in more detail in section 2.5.1 approaches the chattering problem from a fundamentally different angle. Instead of designing the control in continuous time and dealing with discretization effects in a second step, a discrete time description of the design system is adopted.

In this chapter we present a brief study of the sensitivity of super twisting control w.r.t. discrete sampling and first order input dynamics. We then address the design of a constrained integral sliding surface that enables collision avoidance guarantees under bounded control errors. We then present a novel predictive sliding mode control approach extending DSMC that approximates the continuous-time performance of classic sliding mode control while exhibiting heavily improved chattering properties in the presence of discrete sampling and input dynamics.

Super twisting guidance laws with discrete sampling and inner loop dynamics

Encouraged by the ability of super twisting control to provide continuous control without requiring the derivative of the sliding variable, the authors of [START_REF] Galzi | Formation flight strategy and control using higher order sliding mode[END_REF][START_REF] Yamasaki | Sliding mode based integrated guidance and autopilot for chasing UAV with the concept of time-scaled dynamic inversion[END_REF][START_REF] Galzi | Unmanned rotorcraft tight formation flight control using sliding mode control driven by sliding mode disturbance observers[END_REF][START_REF] Yamasaki | Separate-channel integrated guidance and autopilot for a path-following UAV via high-order sliding modes[END_REF][START_REF] Galzi | Closed-coupled formation flight control using quasicontinuous high-order sliding-mode[END_REF][START_REF] Rao | Sliding mode control-based autopilots for leaderless consensus of unmanned aerial vehicles[END_REF] propose the use of STSMC control for aircraft predecessor tracking. Consider the super twisting controller proposed in [START_REF] Galzi | UAV formations control using high order sliding modes[END_REF] for tight formation flight. Extending trivially from 2D to 3D and assuming decoupled load factor tracking loops [START_REF] Bolting | Discrete sliding mode control of small UAS in tight formation flight under information constraints[END_REF], it is given by

ng p = α p |σ p | 1/2 sign(σ p ) + β p sign(σ p )dt (2.5.5)
where p = 1...3 indicates the axes of the guidance frame. We consider the simple LTI sliding surface defined by

σ(t) = Gx(t) (2.5.6) x(t) = ∆p g (t) ∆v g (t) (2.5.7)
where G ∈ R 3×6 . In the following, for the continuous-time case, the dependence on time is dropped for notational convenience. With

G = G 1 G 2 (2.5.8)
the position error dynamics in sliding mode (σ = 0) are

0 = G 1 G 2 ∆p g ∆v g (2.
5.9)

∆v g = -G -1 2 G 1 ∆p g (2.5.10)
Selecting -G -1 2 G 1 as Hurwitz ensures that ∆p asymptotically converges to zero while in sliding mode. From equation (2.4.21) we have the open loop σ dynamics of relative vector degree

r = (1 1 1) T σ = G ẋ (2.5.11) = G ∆v g i + ξ v g ng i+1 + ξ a (2.5.12)
The super twisting algorithm requires bounds on the disturbance magnitude and rate.

To simplify defining these bounds in this purely academic setting, we select an inertial guidance frame aligned with the local inertial frame, making R ge (t) = I 3 and removing the unmatched disturbance ξ v , a constant relative position command and no atmospheric perturbations, reducing the remaining disturbance term ξ a to

ξ a = ve i -g g + gn g w,i+1
(2.5.13)

Pre-compensating known terms by the equivalent control

ng i+1,eq = g -1 (g g -G -1 2 G 1 ∆v g i ) (2.5.14)
reduces equation (2.5.12) to

σ = G 2 g ng i+1,1 -ξa (2.5.15)
where ng i,1 is the SMC control compensating the unknown disturbances ξa = ve i + gn g w,i+1

(2.5.16)

We can easily enforce bounds on ve i by tracking a virtual predecessor following a sinusoidal trajectory satisfying ve i,p ≤ 1 for p = 1, 2, 3. In the absence of atmospheric perturbations, load factor tracking errors captured by the second term are solely due to unknown input dynamics. Assuming that n g i+1 is the output of a stable LTI system of order one, representing well behaved inner loop load factor tracking dynamics

ṅg i+1 (t) = τ -1 (n g i+1 (t) -n g i+1 (t)) (2.5.17)
introduces the obvious issue that equation (2.5.17) breaks the matching condition on top of which classic SMC stability and performance proofs are built. When trying to model this dynamics as feedthrough with additive disturbances, consider the system in the Laplace domain

n g i+1 (s) = ng i+1 (s) 1 + τ s (2.5.18)
The additive disturbance is

n g w,i+1 (s) = 1 - ng i+1 (s) 1 + τ s (2.5.19)
It follows immediately that for high input frequencies (s = jω → ∞), the disturbance n g w,i+1 (s) approaches the commanded input ng i+1 (s), reducing the input margin available for disturbance rejection to zero. Existing proposed HOSM tight formation flight laws sometimes circumvent this fundamental issue by neglecting the rotatory UAS dynamics altogether (e.g. [START_REF] Galzi | Formation flight strategy and control using higher order sliding mode[END_REF]).

Illustration: Standard Vehicle Model closed loop response

To demonstrate the effects of neglecting inner loop dynamics and discrete implementation in the control design model, we apply 2.5.5 to a pair of two SVM models. Control gains are computed according to [START_REF] Shtessel | Sliding mode control and observation[END_REF]. For implementation, the STCSMC is sampled with a zero-order hold scheme, leading to the Time-sampled Super Twisting (TSST) controller . We evaluate the two performance indices -mean squared error and maximum absolute error -over the cruise flight flight segment of the FFB benchmark mission. Note that since only the vertical channel is considered, the corresponding weighting coefficient of the weighted mean squared error (see equation ( 2

Discontinuity in inner loops

A further, more fundamental limitation arises from the discontinuous first derivative of STSMC control that acts as tracking reference for the inner loops. Consider a load factor tracking law (without direct lift control) with the pitching moment coefficient

C mcg = C m 0 + C mα α + C mq q + C m δe δ e (2.5.20)
and assume a hypothetical equivalent control enforcing sliding mode pitch angle command tracking and an LTI first order sliding surface. For clarity of exposition, level flight in calm air and thus θ = q and α = θ is assumed. The sliding surface with time constant τ θ is given by

q = -τ -1 θ (θ -θ c ) (2.5.21) σ = -τ -1 θ (θ -θ c ) -q (2.5.22)
Computing the equivalent control from σ = 0, one obtains:

0 = -τ -1 θ (q -θc ) -q (2.5.23) = -τ -1 θ (q -θc ) - qSc J yy (C m 0 + C mα θ + c 2V a C mq q + C m δe δ e ) (2.5.24
)

δ e = C m δe qSc J yy -1 -τ -1 θ ( θc -q) + qSc J yy (C m 0 + C mα θ + c 2V a C mq q) (2.5.25)
Note that the equivalent control is not only inversely proportional to the desired time constant, it is also proportional to the commanded pitch rates θc . Thus whether magnitude constraints of the elevator actuator are violated depends on the pitch rate reference trajectory, which can be considered as a disturbance in equation (2.5.25). Further dif-ferentiating, we obtain

δe = C m δe qSc J yy -1 -τ -1 θ ( θc -q) + qSc J yy (C m 0 + C mα θ + c 2V a C mq q) (2.5.26)
Note that to guarantee bounded actuator input rates δe , the second derivative of the commanded pitch angle, θc needs to be bounded. Applying super twisting sliding mode control, the first derivative of the control input is discontinuous, making the second derivative unbounded. These unbounded signals that appear as virtual inputs propagate to actual actuator commands constitute a further, more fundamental limitation of applying STSMC at the guidance level.

In the sequel we study the application of discrete sliding mode approaches to extend the performance of TSST control at high sampling rates and fast inner loops towards a larger envelope that covers regions compatible with the avionics and inner loop dynamics of more realistic fixed-wing UAS models such as the FFB Cularis UAS (see section 1.2).

Predictive discrete sliding mode control

The detrimental effects of discrete sampling have given rise to discrete variants of continuous time sliding mode techniques.

Discrete Sliding Mode Control

Discrete Sliding Mode Control (DSMC) techniques such as proposed in [START_REF] Bartoszewicz | Discrete-time quasi-sliding-mode control strategies[END_REF] approach the chattering issue of continuous time SMC by employing discrete-time Lyapunov functions (DTLF) and a discrete-time control design model. The condition on negative definiteness of the time derivative of a continuous time Lyapunov function

V (σ) < 0 (2.5.27)
is replaced by a condition on the DTLF decreasing over two successive samples

V (σ k+1 ) < V (σ k ) (2.

5.28)

A desired convergence rate can be enforced by requiring

V (σ k+1 ) = ΨV (σ k ) (2.5.29)
with a Ψ ∈ (0, 1). The choice of Ψ allows to trade off control effort against reaching time and boundary layer size. DSMC control of this type has been mostly applied to electromechanical systems. See for instance experimental results reported in [START_REF] Chuang | Applying discrete dynamic integral sliding surface control to hydraulic position control[END_REF], confirming the expected chattering alleviation in a hydraulic piston position control appli-cation. Intuitively, the advantages of DSMC when it comes to chattering are due to its being based on a control objective that is closer to our actual objective, which is confining the system states to the sliding surface. Instead of requiring that the derivative of the sliding variable points towards the sliding surface, we require that the state gets closer to it at each sampling instant. This predictive nature allows us to avoid chattering, which is nothing else than a local overshoot over the sliding surface. There are strong parallels between Model Predictive Control and DSMC control and as we will see a DSMC control of the type given by equation (2.5.29) can be cast as a special kind of MPC problem. To the best of the author's knowledge, discrete sliding mode control techniques have not been applied to the problem of aircraft formation flight so far.

The promising properties of DSMC come with a price. To perform the one-step-ahead prediction, a model is required that at least locally well approximates the dynamics of the sliding variable. That being said, as a particular disadvantage, the input gain has to be known, while the continuous time super twisting algorithm only requires its sign and bounds. This property precludes the use of DSMC with inner loops interfaces based on airframe orientation or angle of attack. The particular choice of a load factor interface we make in this work however circumvents this issue, opening the way to apply DSMC techniques.

While heavily mitigating chattering due to discrete sampling, inner loop dynamics negatively affect classic DSMC control by introducing uncertainty into the prediction implicit in the reaching law (2.5.29). In this work we present a extended discrete SMC control scheme that mitigates this chattering effect of inner loop dynamics. We find that established DSMC control laws can be written as equivalent model predictive control problem. This allows us to extend existing algorithms to integrate inner loops states in the prediction step that is solved by established MPC techniques. This also enables us to readily take into account on-line magnitude and rate constraints of inner loops, a further benefit over existing techniques.

In the following, we first derive the discretized predecessor-follower dynamics. We then present an optimization-based procedure for the design of discrete-time sliding surfaces consistent with the constraints of tight formation flight. We then present a predictive discrete-time sliding mode (PDSCM) guidance law and a discrete-time disturbance estimation scheme. Finally, simulations on the FFB benchmark show improved position tracking performance compared to LQ baseline guidance laws.

Discrete predecessor-follower dynamics

The matching condition in discrete time

In real world physical systems, disturbances are always unmatched, due to the inevitable presence of unknown, possibly very fast high order input dynamics. As one example, the usually neglected mechanical link between a servo actuator and its corresponding control surface represents a fairly complex, very fast and stable nonlinear system in itself. Only by the specific choice of the design model the matching condition can be satisfied. Forward Euler discretization (as selected for instance in an earlier version [START_REF] Bolting | Discrete sliding mode control of small UAS in tight formation flight under information constraints[END_REF] of the control presented here, not considering inner loop dynamics) preserves the matching conditions of the continuous time dynamics. Under exact discretization however, inputs that are separated from the disturbances by fast stable input dynamics can become matched due to the decay of these parasitic dynamics within one sampling period. Considering the continuous time dynamics of a pair of UAS as given by equation (2.4.21) and assuming that the inner loop tracking behavior is reasonably well approximated by the first order LTI system given by equation (2.5.17) we have

∆ ṗg i = ∆v g i + ξ v ∆ vg i = gn g i+1 + ξ a ṅg i+1 = τ (n g i+1 -n g i+1 )
(2.5.30)

For greater clarity and ease of notation, we use in the following a generic state space representation of 2.5.30:

ẋ = Āx + Bu Bd u d (2.5.31)
Discretizing 2.5.31, some kind of assumption has to be made about the behavior of control and disturbance inputs between sampling instants. The zero-order hold assumption, most common in MPC algorithms, assumes that inputs keep a constant value between samples. While this can be enforced for control inputs, disturbances generally evolve in between samples and can furthermore only be estimated with some uncertainty. To capture this additional degree of freedom, we split disturbances into a known (correctly estimated) part, constant over the next sampling period, and a bounded but unknown part

d = d + ∆d (
0 3 I 3 0 3 0 3 0 3 gI 3 0 3 0 3 -τ I 3   
(2.5.37)

Bu =    0 3 0 3 τ I 3    (2.5.38) Bd =    I 3 0 3 0 3 I 3 0 3 0 3    (2.5.39)
Note that both disturbances in 2.5.33 are unmatched. Provided the zero-order hold assumption holds, the exact discrete counterpart of 2.5.31 is given by

x k+1 = e ĀTs x k + e ĀTs Ts 0 e -Āt Bdt    u k dk ∆d k    (2.5.40) = Ax k + B    u k dk ∆d k    (2.5.41)
where A, B are the discrete state transition matrix and input matrix respectively and T s is the sampling time.

Constrained discrete integral sliding surface

In quasi sliding mode, the sliding surface with its boundary layer constrains the dynamics of the system states. During the reaching phase, apart from a stability condition that ensures contraction of the sliding variable, no guarantees concerning the evolution of the system state are available. More particularly, additional effort would be required to obtain bounds on the state space region of the system stays within during the reaching phase, i.e. to ensure collision avoidance or envelope protection (e.g. V a > V a min ).

Integral sliding mode approaches avoid the reaching phase altogether by designing a dynamic sliding surface that starts at the initial state. This dynamic sliding surface can be designed by applying any standard control technique such as LQR, H ∞ , pole placement etc. to a nominal design model. In formation flight, state constraints such as collision avoidance regions and minimum and maximum airspeed need to be respected not only during the reaching phase but over the whole mission time. An example that allows to deal with velocity constraints is the nonlinear sliding surface proposed in [START_REF] Shah | Guidance of air vehicles: A sliding mode approach[END_REF]. However, the equivalent control necessary to evolve along the sliding surface additionally needs to allow for a sufficient margin to reject future disturbances unknown at the time of surface computation.

Very recently, combinations of MPC and sliding mode control have been proposed that accommodate this set of requirements appearing in many control problems beyond aircraft guidance in a very natural manner. In [START_REF] Liao | Robust model predictive control with discrete-time integral sliding surface[END_REF][START_REF] Steinberger | Model predictive output integral sliding mode control[END_REF] an MPC control law is augmented with a first order switching SMC term and good performance in the presence of matched disturbances on a test bench is demonstrated while showing the usual input chattering. Here we apply this idea to the problem of tight fixed-wing formation flight, replacing the switched input of classic SMC by a predictive discrete sliding mode input to reduce chattering.

Problem statement

Starting from an initial relative state x 0 , we ask for a sliding surface (SSF) that satisfies σ(x 0 ) = 0 (2.5.42)

i.e. there is no reaching phase, and u eq,k ∈ U eq ∀k ≥ 0 (2.5.43)

where u eq,k is the equivalent control, i.e. the commanded load factors required to have the nominal, disturbance-free system follow the SSF, U eq is an appropriately tightened load factor constraint set that leaves sufficient margins to reject future disturbances while following the sliding surface. Note that the sequence u eq,0 ...u eq,N -1 directly provides the equivalent control necessary to follow the SSF in the nominal disturbance free case. This constraint tightening approach is known from robust MPC, see e.g. [START_REF] Mayne | Robust model predictive control of constrained linear systems with bounded disturbances[END_REF]. Furthermore we require

x * k ∈ X ∀k ≥ 0 (2.5.44) x * N ∈ X N (2.5.45)
i.e. the state trajectory x * k representing the integral SSF is required to respect constraints embodied by the state constraint set X e.g. the velocity must not undercut a given value related to the stall speed and the position must not coincide with the space occupied by other members of the formation. The state constraint set can furthermore be tightened, taking into account maximum control errors of the sliding mode control that drives the system along the SSF, to guarantee that the original constraint set is not violated under disturbed feedback control. As an example, consider figure 2.6 where the SSF respects constraints on vertical velocity, while a control error of about 0.3m/s causes constraint violation in flight in turbulent air. By tightening state constraints, robust constraint compliance can be enforced, see figure 2.7.

Furthermore, the state is supposed to be driven into a neighborhood of the sweet spot in the predecessor's wake, represented by the terminal state constraint set X N . The size of X N depends on the region of attraction of the peak seeking algorithm employed to converge to the position of maximum upwash gain and the amount of uncertainty w.r.t to it.

The problem can be posed as a standard MPC problem with terminal constraint set and prediction horizon N * min u eq,0 ...u eq,N * -1 and solved as a quadratic program using fast embeddable QP solvers such as [START_REF] Mattingley | Cvxgen: A code generator for embedded convex optimization[END_REF].

x * 1 ...x * N * N * -1 k=0 u T eq,k Ru eq,k + x * T k+1 Qx * k+1 subject to u eq,k ∈ U eq x * k ∈ X x * N * ∈ X N x * k+1 = Ax * k + B u u eq,k ∀k = 1...N * (
Horizon length It is desirable to minimize the prediction horizon length for computational reasons. Secondly, while moving along the SSF to the sweet spot, no energy gains can be enjoyed by the follower, adding an additional motivation to find a minimumtime solution. One answer to the minimum time horizon length is given by the notion of reachable sets. If the reachable set corresponding to N contains X N , X N is reachable for an horizon of N. Computing reachable sets is computationally involved even for LTI systems. The minimum horizon length has however to be determined on line, since X is partly unknown and time-varying due to a-priori unknown positions of other formation members.

In their recent contribution [START_REF] Van Den Broeck | Model predictive control for timeoptimal point-to-point motion control[END_REF], the authors propose finding the minimum horizon N * for time-optimal model predictive motion control by sequential feasibility checks of the corresponding quadratic program. An implementation example achieving solution update rates of 200 Hz on a standard PC is given, suggesting that this approach is a promising candidate for the minimum time constrained integral sliding surface (CISSF) problem, given the powerful computational resources of today's UAS. It has been implemented in our simulations to identify N * .

From standard to generalized PDSMC

We proceed now to designing predictive discrete sliding mode laws that have the UAS follow the constrained sliding surface. Employing the CISSF defined in 2.5.3, the open loop sliding variable dynamics follow to

σ k+1 = G(x k+1 -x * k+1 ) (2.5.47) = G(Ax k + B u B d B d    u k dk ∆d k    -Ax * k -B u u eq,k ) (2.5.48)
Applying the equivalent control sequence resulting from solving 2.5.46

σ k+1 = G    A (x k -x * k ) B u B d B d    ūk dk ∆d k       (2.5.49)
where

ūk = u k -u eq,k ∈ U ∀k ≥ 0 (2.5.50)
is the control available to the sliding mode control for disturbance rejection.

Existing DSMC algorithms such as [START_REF] Bartoszewicz | Discrete-time quasi-sliding-mode control strategies[END_REF] seek to enforce a prescribed decay rate of σ. Due to the unknown disturbance component d, ideal sliding mode is not achievable for finite sampling rates. It is however possible to drive the system into quasi-sliding mode, defined by

|σ k | ≤ ∀k ≥ 0 (2.5.51)
and it can be shown that the boundary layer thickness depends quadratically on the sampling time [START_REF] Bolting | Discrete sliding mode control of small UAS in tight formation flight under information constraints[END_REF][START_REF] Su | An O(T 2 ) boundary layer in sliding mode for sampled-data systems[END_REF] and linearly on the prescribed decay rate. Thus in order to minimize the boundary layer thickness, a small Ψ is desirable. At the same time, for static sliding surfaces, the initial system state may be far off the sliding surface, limiting the choice of Ψ to avoid control input saturations.

PDSMC

To meet these conflicting requirements, a time-varying, less restrictive reaching law is proposed in this work, leading to what we denote Predictive Discrete Sliding Mode Control (PDSMC). It is conceptually similar to what the authors of [START_REF] Houda | New discrete predictive sliding mode control for non minimum phase systems[END_REF][START_REF] Perez | Design of an explicit constrained predictive sliding mode controller[END_REF] proposed in the context of chemical processing. In contrast to these existing approaches though, we take hard constraints on input magnitude and rate saturations into account and focus on very short prediction horizons. As the linear reaching laws of type (2.5.29), we enforce a contraction of σ towards the sliding surface at each time step. In contrast to the reaching law (2.5.29), or the optimization based control proposed in [START_REF] Perez | Design of an explicit constrained predictive sliding mode controller[END_REF] the step towards the sliding surface is maximized by solving at each time step the minimization problem min

ūk V (σ k+1 ) subject to ū ∈ U ⊂ R p ∀k ≥ 0 (2.5.52)
Note that this control law has no free tuning parameters that need to be chosen by the control designer but is defined by hard input constraints dictated by the inner loops . It is equivalent to choosing a smaller reaching matrix Ψ when closer to the sliding surface, leading to tighter bounds on the boundary layer thickness and thus reducing the maximum tracking error. A simple choice for the discrete Lyapunov function

V (σ) is V (σ) = 1 2 σ 2 .
Inner loop input rate saturations can be taken into account by requiring

ūk -ūk-1 ∈ U ∆ (2.5.53)
which can be enforced by setting

Ūk = Ū ∩ (ū k ⊕ U ∆ ) (2.5.54)
where ⊕ is the Minkowski sum of two sets given by A ⊕ B = {a + b|A ∈ A, b ∈ B}, i.e. we limit the magnitude constraint set to its part that is consistent with rate constraints. The optimization problem 2.5.52 can efficiently be solved as a small quadratic program, just like linear MPC problems and the same available fast QP solvers can be employed. PDSMC control shares more similarities with MPC and can be written as a special MPC problem. This property is quite remarkable, taking into account PDSMC originally being based on a discrete version of sliding mode control. As a convenient benefit, established ideas and techniques from the field of MPC can readily be exploited in the design of PDSMC control, such as stability constraints.

Stability constraints

Terminal state constraint sets are the dominant way to encode stability into MPC problems, (where usually N > 1) either as dual mode MPC or based on contraction. In dual mode MPC the existence of an infinite-horizon unconstrained controller is assumed, e.g. an LQR whose region of attraction (RoA) can be computed. A terminal state constraint set that is an exact representation or outer approximation of this RoA ensures that the last state of the prediction horizon lies within this set. At the end of the prediction horizon the unconstrained controller can thus take over and drive the system to the origin. In contraction based MPC the terminal state constraint set enforces that a Lyapunov function of the state has decreased. Contraction based stability constraints generally allow for shorter prediction horizons especially when the system is farther away from the origin. The same principle could be applied to PDSMC by including a linear equivalent of the contraction constraint

σ 2 k+1 < σ 2 k (2.5.55) as σ k+1    ∈ [0, σ k ), if σ k ≥ 0 ∈ (σ k , 0], if σ k < 0 (2.5.56)
This constraint can be used to check for the feasibility of the optimization problem (OP) 2.5.52. However, if feasible, the constraint 2.5.55 is implicitly satisfied by solving 2.5.52 and does not need to be implicitly included into the QP, reducing its complexity and speeding up solution times.

It is worthwhile at this point to take a step back and put this control as well as classic DSMC approaches of type (2.5.29) into the context of model predictive control to appreciate its particular properties. From an MPC perspective, the PDSMC control solves a state trajectory tracking problem with an extremely short prediction horizon of N = 1 (simulations suggest that sampling times of the prediction model of the order of the inner loop time constant provide good performance). This is a particular case of contraction-based MPC, brought forward e.g. in [START_REF] De Oliveira Kothare | Contractive model predictive control for constrained nonlinear systems[END_REF][START_REF] Alamir | Contraction-based nonlinear model predictive control formulation without stability-related terminal constraints[END_REF]. Contraction based stability in MPC typically allows for much lower prediction horizons than other classic approaches to ensure stability of MPC algorithms, such as dual mode MPC However, minimum prediction horizons for stability reported in the MPC literature typically exceed one.

It is insightful to consider that the controlled system, in spite of the discrete-time prediction model, evolves in continuous time. A minimum horizon required for stability by contraction exist therefore as a minimum prediction time horizon. By choosing a sufficiently large sampling time exceeding this minimum prediction time horizon, contraction over one sampling time step becomes possible. This is an assumption implicitly made in existing DSMC approaches such as [START_REF] Bartoszewicz | Discrete-time quasi-sliding-mode control strategies[END_REF]. Considering a long sampling time under zero-order input hold is equivalent to an extreme case of what is denoted as "move blocking" in MPC. That being said, posing the PDSMC control problem in an optimization framework also allows for future interesting modifications such as prediction horizons N > 1 in conjunction with shorter sampling times. The OP 2.5.52 minimizes the maximum absolute value of the sliding variable at each sampling step, in effect minimizing the boundary layer. By computing an upper bound on σ 2 k+1 , the maximum boundary layer thickness can readily be obtained. This is an optimal, i.e. minimum achievable boundary layer thickness under the given input constraints and not a design parameter, in contrast to SMC approaches based on smoothing the CSMC discontinuity by e.g. saturation or sigmoid functions. This sets PDSMC control apart from existing DSMC schemes where the boundary layer is function of the decay rate and as such needs to be chosen by the control designer. Extending the prediction horizon further than one, the question arises of how to maintain this property. The LQ type costs used in most MPC algorithms form a weighted sum of control errors and as such are not suited. H ∞ MPC control approaches based on Linear Matrix Inequalities have been proposed [START_REF] Orukpe | Model predictive control based on mixed h2/h control approach[END_REF] and adapting this strategy for PDSMC control seems the natural way to approach this problem. However, while mature and fast embedded QP solvers based on automatic C code generation are readily available, this appears not to be the case for LMI solvers [START_REF] Vandenberghe | Algorithms and software tools for LMI problems in control: An overview[END_REF].

We present in the following a different approach based on contraction constraints that can be solved as a standard QP.

Generalized PDSMC

Consider the following optimization problem

P σ with σ k = Gx k min ū0 ...ū N -1 x 1 ...x N x T 1 (G T PG)x 1 + N -1 k=0 ūT k Rū k subject to ūk ∈ U Gx 1    ∈ [0, Gx 0 ), if Gx 0 ≥ 0 ∈ (Gx 0 , 0], if Gx 0 < 0 Gx k    ∈ [0, Gx 1 ], if Gx 0 ≥ 0 ∈ [Gx 1 , 0], if Gx 0 < 0 ∀k = 2...N x k+1 = Ax k + B u ūk + B d d0 ∀k = 0...N -1
Assuming feasibility, P σ enforces contraction of the first predicted sliding variable σ 1 . It enforces furthermore that the state trajectory does not cross the sliding surface and that all following σ k do not exceed σ 1 , thus in effect minimizing the predicted boundary layer. The novelty of this approach consists in optimizing over the next predicted state, and using this predicted state to constrain the following ones. The matrix R provides, as in standard MPC problems, a degree of freedom to the control designer to trade off control effort against boundary layer thickness.

Note that the previously considered DSMC reaching law can be stated in the same optimization framework by replacing contraction inequality constraints by equality con-straints [START_REF] Perez | Design of an explicit constrained predictive sliding mode controller[END_REF]. Various formulations are possible, one example is min

u k ,x k+1 ūT k ūk subject to x k+1 = Ax k + B u ūk + B d dk Gx k+1 = ΨGx k (2.5.57)
In the following we direct our attention towards the minimum horizon single input PDSMC control 2.5.52 which is a special case of P σ with N = 1, p = 1 and R = 0. Due to the very small problem size it can actually easily be solved analytically by what amounts to a minimal version of an active set QP algorithm. Considering

σ k+1 = G(Ax k + Bu k + B d dk + B d ∆d k ) (2.5.58) V = 1 2 σ 2 k+1
(2.5.59)

At an unconstrained minimum we have δV δu = 0, thus for the nominal control (∆d = 0)

δV δu = (G(Ax k B d dk ) + GBu k )GB (2.5.60)
= 0 (2.5.61)

and the unconstrained PDSMC control is

u * = -(GB) -1 (GAx k + GB d dk ) (2.5.62)
If the unconstrained control u * ∈ U , it is identical to the constrained control. Otherwise, we compute V at the boundaries of U and select the smallest one, e.g. for a box constraint U = [u, ū], we obtain two values of V (σ) corresponding to applying u or ū. Note that, interestingly enough, for a linear sliding surface of type σ = Gx and a linear system, in a neighborhood of the sliding surface where input constraint are inactive, the single input PDSMC control with N = 1 results in a linear state feedback plus a disturbance rejection term. The non-iterative computation of u * enables fixed-time implementation, an interesting general property for flight control, where computing time guarantees need to be given. We have not commented so far on how to obtain the disturbance estimate dk . In the following we investigate suitable real-time estimators.

Disturbance estimation

Prediction based PDSMC control builds on the assumption that some kind of estimate of unknown disturbances due to predecessor motion, the rotating guidance frame etc. is available. In this section we present a continuous time sliding mode observer and a discrete disturbance observer and weigh their respective merits.

=gn g i+1 + R ge (t) ve i -Ω g (t)(∆v g i + ∆ ṗg c,i (t)) -Ω g (t)∆ ṗg -Ωg (t)∆p g -∆p g c,i (t) + g g + gn g w,i+1
(2.5.66)

Disturbance elements such as those due to gravity and imperfect load factor tracking acting on each vehicle due to atmospheric turbulence and predecessor wake vortices can in principle be estimated using on board accelerometers. Disturbances due to predecessor motion, require high-rate communication between vehicles.

Super twisting sliding mode disturbance observer

A different approach can be taken by estimating the sum of disturbances with a sliding mode observer (SMO). The problem tackled here is basically a problem of estimating unknown acceleration inputs by real time signal differentiation, in the presence of discrete sampling and measurement noise. Sliding mode observers, just like sliding mode controllers, have very attractive properties in a continuous time framework. SMO techniques have been widely applied to disturbance estimation problems in the context of fault detection and real time signal differentiation. The basic idea is to treat the estimation problem as a control problem and to drive the error between the velocity of a model of a pair of vehicles internal to the observer and the estimated actual relative velocity to zero by means of a virtual control input. Once the error has converged to zero, the virtual control input generated by the SMO is an estimate of the accelerations acting on the system. Consider the vertical channel of 2.5.66

∆ vg i,z = a z (2.5.67)
where ∆ vg i,z is the relative velocity of the observer model, and a z is the virtual control input. Defining the sliding variable σ and its first time derivative σ = ∆v g i,z -∆v g i,z

(2.5.68)

σ = gn g i+1,z + ξ a -a z (2.5.69)
we see that σ is of relative degree one and both classic first order SMC or second or-der super twisting control can be applied. First order SMC generates high-frequency switching in the virtual control ñ and low-pass filtering needs to be applied to obtain a smooth estimate, introducing additional lag in the observer response. On the other hand, the super-twisting algorithm readily generates continuous virtual inputs, it is therefore the natural choice

a z = α|σ| 1/2 sign(σ) + β sign(σ)dt (2.5.70)
where α, β have to be tuned according to bounds on ξ a , ξa , e.g. to satisfy conditions given in [START_REF] Shtessel | Sliding mode control and observation[END_REF]. Applying this observer to the relative velocity between two UAS in the guidance frame, we obtain an estimate of its total time derivative. To recover the unknown disturbance, consider equation (2.5.66)

ξasmo,k = a z,k-1 -gn g i+1,k-1
(2.5.71)

Discrete disturbance observer

The finite difference differential is given by

∆ vg i,z,k = ∆v g i,z,k -∆v g i,z,k-1 T s (2.5.72)
and an estimate of the unknown disturbance is recovered by the discrete disturbance observer (DDO) as

ξa ddo ,k = ∆ vg i,z,k -gn g i+1,k-1 (2.5.73)
The sensitivity of 2.5.72 w.r.t. additive noise w k as function of sample time can easily be analyzed by considering simulation, a synthetic disturbance ξ a is injected. This way magnitude and derivative bounds are perfectly known and the gains of the SMO can be computed without any assumptions. Measurement noise corresponding to estimation quality achievable with GNSS carrier phase time differencing [START_REF] Freda | Time-differenced carrier phases technique for precise GNSS velocity estimation[END_REF] is added to the velocity signal. See [START_REF] Smaoui | A robust differentiator-controller design for an electropneumatic system[END_REF] for a recent application of SMO for the piston position control of an electropneumatic system. The choice of SMO over finite-difference differentiation is backed in [START_REF] Smaoui | A robust differentiator-controller design for an electropneumatic system[END_REF] by a simulation example showing smaller errors in the estimated signal derivative in the presence of noise. Our simulations confirm this superiority for fast sampling. Taking differentiator sampling time into account, the picture becomes more involved. While at high sampling rates the SMO handles even noisy velocity observations well, derivation performance deteriorates with lower sampling rates, similar to sampled continuous time sliding mode control performance. See figures 2.8, 2.9 for snapshots of observation performance at sampling rates of 10 -2 s and 10 -1 s respectively. In the first case, closer to its continuous time design conditions, the SMO shows significantly smaller estimation errors. For a slower sampling rate of 10 -1 s, the picture reverses, the disturbance estimate of the SMO is dominated by numerical chattering, while the DDO estimate follows well the true disturbance, with the phase lag inherent to finite-difference differentiation, see figure 2.9. This small case study illustrates how sliding mode observer performance, just like its SMC dual, is tied to its discrete implementation and can be outperformed by very simple schemes under unfavorable conditions. See figure 2.10 for a sweep of the rms and maximum absolute error of the SMO and the DDO over a range of sampling rates. The DDO differentiation error in the noisy case evolve as indicated by equation (2.5.75) and converges to the noise-free error for lower sampling rates. For low sampling rates, the error due to phase lag dominates, and the error in case of noisy signal converges to that corresponding to the noise-free case. The sliding mode observer shows superior performance only for fast sampling rates down to about 50 Hz.

∆ vg i,z,k + e = ∆v g i,z,k + w k -∆v g i,z,k-1 -w k-1 T s (2.5.74) = ∆ vg i,z,k + w k -w k-1 T s ( 2

A preliminary closed-loop evaluation

We apply PDSMC control and the discrete disturbance observer 2.5.

FFB simulations

We first consider one individual UAS tracking a virtual leader on the benchmark trajectory. Longitudinal and lateral spacing is maintained by baseline LQ control. The trajectory speed in all scenarios is 15 m/s. Note that position tracking errors are normalized with span width and commanded load factors as well as elevator and flap commands with their maximum admissible values. A low controller sampling rate of 10 Hz is selected, compatible with observation rates of low-cost GNSS receivers. Running guidance laws at the same rate as GNSS observations become available, no high-rate exchange of accelerometer observations and no IMU/GNSS data fusion is required, lowering computational load.

Turbulent air

In this scenario, planar turbulent wind of ||v e w,a || = 3 m/s aligned with the cruise flight leg of the FFB trajectory is simulated. Tracking errors of the PDSMC control without disturbance estimation, displayed in figure 2.15) are slightly worse than those achieved by the LQ baseline laws, see figure 2.14). Introducing the disturbance estimation scheme presented in section 2.5.5 appreciably reduces especially peak errors by about 50 % during predecessor maneuvers.

Large initial errors

Large initial position errors are readily accommodated by the predictive integral sliding surface, see figure 2.17. Peak load factors are of similar magnitude as in the case of small

Calm air

Flying the benchmark in calm air confirms the importance of disturbance estimation to suppress tracking errors due to predecessor motion.

Scalability

When adding more members to the formation, linear control with local state feedback is known to lead to position errors growing with vehicle index, i.e. mesh instability [START_REF] Seiler | String instabilities in formation flight: Limitations due to integral constraints[END_REF]. As shown in chapter 2, the PDSMC for N = 1 boils down to linear control, leading us to expect mesh instability when applying it to larger formations. This expectation is confirmed by simulations of three vehicles, see figure 2 Chapter 3

Set membership localization

The position of the sweet spot in the wake of a preceding UAS allowing for maximum range enhancement is airframe-dependent and only approximate models of the complex flow field forming the wake are available. Peak-seeking algorithms, see e.g. [START_REF] Brodecki | Autonomous formation flight control system using in-flight sweet-spot estimation[END_REF], have been proposed to deal with this uncertainty. Peak seeking algorithms attempt to minimize on-line some cost function related to the wake energy gain to converge to the point of minimum power consumption without precise knowledge of its position. Note however that their gradient based nature relies on convexity of the cost, which can exhibit local minima, see section 1.2. That being said, some exploration of the wake might be necessary to identify a suitable starting point for peak seeking that has a high probability of leading to convergence to the global maximum of wake energy gain. During this exploratory phase, while peak seeking is active and while flying in a static tight formation, collision avoidance needs to be ensured under localization uncertainty as well as bounded control errors.

Accuracy and safety

Localization errors generally have a greater detrimental effect on formations of small UAS than larger ones, since admissible guidance errors scale with airframe size. Using COTS consumer grade GNSS hardware, relative localization errors in the centimeter range are possible using readily available Real Time Kinematic (RTK) algorithms, see e.g. [START_REF] Gross | Robust UAV relative navigation with DGPS, INS, and peer-to-peer radio ranging[END_REF], or commercial low cost RTK systems such as [START_REF] Emlid | Emlid Reach RTK[END_REF], corresponding to the high precision of GNSS carrier phase observations. This accuracy falls well within the bounds required for TFF even for sub-meter class UAS. Alternative GNSS-independent approaches such as machine vision have been considered but generally provide less accuracy, among other reasons owing to attitude estimation errors of low-cost AHRS, affecting the critical transformation of observations in the camera frame to the local inertial frame, see e.g. [START_REF] Wilson | Vision-aided guidance and navigation for close formation flight[END_REF].

Availability is a major concern when it comes to using GNSS for guidance, and can be prohibitive in close-to ground operations due to multipath and signal obstruction. However, in energy-efficient planar formation flight, excellent availability conditions can be expected. For civil UAS applications, where furthermore no intentional GNSS jamming or spoofing is to be expected, localization based on single frequency RTK systems appears to be well suited, and localization accuracy is thus not a fundamental issue when it comes to tight formation flight of small UAS. However, to ensure collision avoidance, the unavoidable error of relative position estimates needs to be taken into account.

Problem statement

To this end, not only point estimates of relative position vectors are required, but rather an outer estimate of the regions that can be guaranteed to contain the other members of the formation. We want furthermore this estimate to be useful, i.e. only as large as necessary and relying only on sensors realistically available on board a typical small low cost UAS. What is more, the employed algorithm needs to be compatible with the still somewhat limited computational resources of today's small UAS.

Assumptions and constraints

Hardware resources We assume that in view of cost and mass constraints affecting in particular small UAS for civil applications, the set of observations and state estimates available for localization is limited to

• The image stream of a single camera, fixed in the body frame

• Attitude estimates and observations of angular rates, accelerations and the earth's magnetic field in the body frame provided by a MEMS based AHRS

• Single-frequency GNSS carrier phase, Doppler and pseudo range observations

• Ultra-Wide Band (UWB) inter-UAS ranging observations

In this chapter we shortly review existing relative localization approaches. We then present a modified Extended Set Membership Filter (ESMF) compatible with the assumed limited set sensor hardware. Care has been taken to replace iterative operations such as numerical optimization by closed-form approximations to enable fixed-time implementation.

We compare the filter to a standard Extended Kalman Filter to illustrate the benefits arising from the deterministic nature of set membership filtering.

Airborne relative localization: state of the art

Existing relative localization approaches rely mostly on Real-Time-Kinematics (RTK) GNSS positioning, machine vision or a combination of both.

Differential GNSS Using differential GNSS observations, Real Time Kinematics (RTK) systems provide, after ambiguity resolution, estimates of the position vector between two GNSS receivers with the centimeter-level accuracy of carrier phase positioning. GNSS RTK systems usually estimate the unknown carrier phase ambiguities employing pseudo range double differences in what is denoted the "float solution". The estimator of choice is usually a Kalman filter for its ability to estimate the covariance of its estimation error. Once the covariance of the float solution falls within a predefined threshold, a search algorithm is employed that finds with high probability the correct set of carrier phase ambiguities. In that case, the position error is reduced to centimeter level. Two properties of this standard approach are problematic when it comes to flight safety. First, before ambiguities are fixed, only probabilistic bounds on the estimation error of the float solution are available. Second, once ambiguities are fixed, a small nonzero probability that the fixed set of ambiguities is not the actual one remains, due to the statistical nature of the fixing algorithm.

UWB assisted RTK Very cost effective UWB ranging modules based on time of flight observations have emerged on the consumer market over the last years. Their decimeter level accuracy makes them a natural candidate to enhance differential GNSS localization both on the pseudo range level and for RTK solutions. After early experiments for static surveying applications [START_REF] Macgougan | Multiple UWB range assisted GPS RTK in hostile environments[END_REF], the authors of [START_REF] Gross | Robust UAV relative navigation with DGPS, INS, and peer-to-peer radio ranging[END_REF] report robustness enhancements of a conventional RTK solution based on the LAMBDA [START_REF] Teunissen | The LAMBDA method for fast GPS surveying[END_REF] algorithm for relative localization of two UAS in flight. In this approach, UWB range observations are employed in forming the float ambiguity estimate. UWB ranging has as well seen a rising interest for local ground-based localization systems, see e.g. [START_REF] Ledergerber | A robot self-localization system using one-way ultra-wideband communication[END_REF][START_REF] Kempke | Polypoint: Guiding indoor quadrotors with ultra-wideband localization[END_REF].

Machine vision

Both techniques purely based on machine vision [START_REF] Dell'aquila | Real-time machine-vision-based position sensing system for UAV aerial refueling[END_REF] as well as combinations with differential GNSS [START_REF] Mammarella | Machine vision/GPS integration using EKF for the UAV aerial refueling problem[END_REF] have been employed in the community for accurate relative localization of UAS and centimeter level precision is reported in experimental settings, again employing an Extended Kalman Filter. Machine vision approaches suffer from the fundamental limitation that tracked UASs need to be within the field of view of the tracking UAS for each filter update, necessitating multiple of gimbaled cameras in larger formations as well as trajectory planning algorithms that take into account maintaining continuous visual contact. Maintaining visual contact to multiple neighbors becomes more difficult with decreasing separation, as in tight formation flight close-by UAS might obstruct the view towards other ones.

Algorithms Relative positioning typically involves nonlinear observation equations. This is the case also for GNSS double difference observations, due to small linearization errors usually assumed to be linear in the separation vector (see section 3.4.1 where we derive upper bounds on the linearization error and consider also the impact of standalone position and ephemeris uncertainty). The Extended Kalman Filter (EKF) and the Unscented Kalman Filter (UKF) are the predominant state estimation algorithms employed. The Kalman Filter can indeed provide an ellipsoidal confidence set that is guaranteed to contain the true state with a selected probability given that its underlying assumptions of normally distributed, zero mean modeling errors and measurement errors are met. However, although the normal distribution in many cases provides a reasonable model (see figure 3.5), as is well known, in reality both error sources can only be approximated by it and filter consistency needs to be ensured by careful tuning. Often no attempt to identify the probabilistic properties of the linearization error is made, and the filter is tuned in a grey-box fashion. What is more, Kalman filters, being based on probabilistic error models, provide estimates of probabilistic state estimation error bounds. Even given a consistent estimate, probabilistic bounds can however only guarantee probabilistic collision avoidance. Guaranteed bounds of estimation errors are however of crucial importance for flight safety. Set membership filtering is concerned with the problem of providing tight outer approximations of these bounds.

Set Membership Filters for State Estimation

Set membership state estimation filters are the deterministic counterpart to probabilistic filtering frameworks such as the various flavors of the Kalman filter and have received considerable interest for applications that require guaranteed bounds on state estimation errors. They deal with the problem of finding mathematically tractable outer approximations of the set of states compatible with all past observations of the state with their respective error bounds and the constraints imposed by a state transition model with bounded prediction error. An exact representation of this set is generally difficult or impossible to compute and can take on very complex shapes over time.

Assume an unknown state x k ∈ R n contained in the guaranteed state set D k . The states in D evolve according to

x k+1 = f (x k , u k , w k ) (3.3.1)
where w k ∈ W k is some bounded prediction error. Using equation (3.3.1), at each discrete sampling instant the predicted guaranteed set D + k can be computed, usually expanding D k-1 . This is analogous to the state estimate covariance prediction step of a Kalman Filter. Assume now that some sensor provides observations

y k ∈ R p y k = h(x k , e k ) (3.3.2)
with a bounded measurement error e k ∈ E k . From y and E an observation set Y can be computed, given by the set-valued observation equation

Y k = H(D Y,k ) (3.3.3)
where the consistent state set is

D Y,k = {x|h(x) ∈ Y k } (3.3.4)
i.e. the set of all the states that yield an observation inside the observation set Y. Whenever new measurements y come in, the intersection of the current predicted state set D + with D Y provides an updated state set, analogous to the Kalman filter measurement update. This is the standard set membership filtering problem.

As exact representations of D k are usually computationally intractable, numerous approaches have emerged that employ various outer set approximations of lower complexity, rendering the set sum and intersection operations that form the basis of set membership filtering actually implementable. Application of set membership filtering to the problem of relative localization of UAS does not appear to have received a lot of attention so far, especially in contrast to localization schemes involving Kalman filtering. Two examples of related problems are given in [START_REF] Garcia | Guaranteed Quadrotor Position Estimation Based on GPS Refreshing Measurements[END_REF], where GPS standalone position measurements are considered, comparing interval and zonotope set representations, and [START_REF] Gu | ESMF based multiple UAVs active cooperative observation method in relative velocity coordinates[END_REF], studying a cooperative ellipsoid set membership filter for joint target observation by a group of UAS. Interval membership localization filtering based on standalone pseudo range observations is presented in [START_REF] Drevelle | High integrity GNSS location zone characterization using interval analysis[END_REF] for ground vehicle localization. An interval representation is considered as well in [START_REF] Caiti | Localization of autonomous underwater vehicles by floating acoustic buoys: a set-membership approach[END_REF] for set membership localization filtering of underwater vehicles. The authors of [START_REF] Bonnifait | Cooperative Localization with Reliable Confidence Domains Between Vehicles Sharing GNSS Pseudoranges Errors with No Base Station[END_REF] present a cooperative set membership localization algorithm for ground vehicles based on shared pseudo range observations. Differential GNSS observations as well as inter-UAS ranging observations generally lead to arbitrarily oriented guaranteed state sets that can be largely overapproximated by intervals. Thus set representations with additional degrees of freedom are desirable to obtain less conservative constraints for trajectory planning, freeing up more maneuvering space. What is more, tight set approximations can be decisive for the success of a mission relying on tight formation flight -if relative position set estimates are too large, the members of the formation possibly cannot reach their assigned energy-efficient position in the wake. Set membership filters employing ellipsoids [START_REF] Scholte | Online nonlinear guaranteed estimation with application to a high performance aircraft[END_REF][START_REF] Scholte | A nonlinear set-membership filter for on-line applications[END_REF] have been proposed, and even more flexibility at the cost of computational complexity can be obtained by bounding the guaranteed state set by zonotopes [START_REF] Alamo | Guaranteed state estimation by zonotopes[END_REF][START_REF] Alamo | A set-membership state estimation algorithm based on DC programming[END_REF], a special class of polyhedra. In this work we present an adaptation of the Extended Set Membership Filter [START_REF] Scholte | A nonlinear set-membership filter for on-line applications[END_REF] based on ellipsoidal set representations for the relative localization problem. The ESMF is particularly interesting in this context for its direct comparability with the widely popular Extended Kalman Filter.

Remark The major benefit of set membership filtering is the absolute nature of the given guarantees. It relies on our assumed ability to guarantee bounds on observation errors and system modeling errors. In reality, these error bounds are often obtained from a finite amount of observation samples. Assuming even a bounded -thus staying in a deterministic framework -non-uniform probability distribution of the error, this finite observation window might be too small to detect large errors of low probability. Thus, in a strict sense no guarantee, but a high confidence can be given that all future samples will respect the bounds thus obtained, interestingly giving a probabilistic aspect to this so decidedly deterministic class of filtering.

In what follows, after providing some essential notation, we derive bounded-error linearizations of differential GNSS code phase and carrier phase observations as well as of inter-UAS range observations. We then present a modified Extended Set Membership filter that exploits these three observations to provide an ellipsoidal outer approxima-tion of the set containing the relative position vector between UAS. We compare the modified ESMF to a standard EKF in simulation.

Definitions and Notation

In this chapter we make use of interval as well as ellipsoid representations of sets. Some relevant notions and notations are given in the following for ease of exposition.

Intervals Apart from regular scalar real valued intervals we encounter interval vectors as well as interval matrices in this chapter. We use a general interval matrix notation

[x] for the family of matrices x ∈ R n×m where x i,j ∈ [x i,j , x i,j ] for i = 1...n, j = 1...m to cover all three cases. For n > 1, m = 1, [x] is an interval column vector, for n = 1, m = 1,

[x] is an interval scalar. Wherever additive errors are assumed, it can be more convenient to write intervals in center-range notation as

[x] = x + [w x ] where x is the center point of [x] and [w x ] = [x] -[x, x] where [x,
x] denotes a degenerate interval containing only the point x.

Wherever regular terms and interval terms appear in the same expression, such as

[c] = a + [b]
it is assumed that the regular terms (a) represent degenerate intervals ([a, a]) and all operations in the expression are interval operations.

Interval extensions

Evaluating the result of a function over interval arguments amounts to constructing an interval extension of this function. Interval extensions of various basic operators and functions (+, -, exp, log etc.) are known to be sharp inclusions of the real result sets. When constructing interval extensions of more complex functions, care needs to be taken, as alternative but equivalent formulations of the same function can lead to different, possibly heavily conservative result intervals. Take as a purely illustrative example the functions

f 1 (x) = x -x (3.3.5) f 2 (x) = x(1 -1) (3.3.6)
Whereas f 1 , f 2 are equivalent, evaluating both over the interval [x] = [-1, 1] leads to vastly different results:

[f 1 ([x])] = [2, 2] (3.3.7) [f 2 ([x])] = [0, 0] (3.3.8)
This fundamental issue is linked to a variable appearing more than once in an expression. When evaluating the function using interval arithmetic, the dependency between both occurrences -the fact that both represent the same variable -cannot be taken into account, resulting in inflated result intervals. There is no systematic way to avoid dependency issues for general functions, except reducing their impact by skillful rearrangement of terms. As a useful exception, for (piecewise) monotonic functions dependency issues can be circumvented altogether. Monotonic functions are either increasing or decreasing over the interval range of their argument and it suffices to evaluate the corner points of an argument interval and to select the results that form the widest interval. For functions with vector argument {y = f (x)|x ∈ R n , y ∈ R}, monotonicity can be verified on a by-element basis [START_REF] Hansen | Global optimization using interval analysis: revised and expanded[END_REF]. If the elements of the gradient of f (x) do not change sign over a given input interval vector range, the function is monotonic w.r.t. x and it suffices to evaluate the 2 n corner points of [x] to find a sharp inclusion of [f (x)]. Monotonicity is exploited in section 3.5 to find sharp inclusions of linearization errors.

Ellipsoids

The ESMF relies on ellipsoidal set approximations, and the probability level sets of the KF covariance matrices form ellipsoids as well. An ellipsoid in R n is given by

Ω(x c , P) = {x|(x -x c ) T P -1 (x -x c ) ≤ 1} (3.3.9)
with x ∈ R n some point in space, x c ∈ R n the ellipsoid center and a positive definite matrix P ∈ R n×n defining half axes and orientation. Useful for analyzing ellipsoid size metrics and displaying the ellipsoid surface in 3D, equation (3.3.9) can be decomposed into a diagonal matrix P 0 ∈ R n×n (corresponding to an axes-aligned ellipsoid) carrying the squares of the half axes on its diagonal, and a matrix R ∈ R n×n performing the rotation to its original attitude

(x -x c ) T (R T P 0 R) -1 (x -x c ) ≤ 1 (3.3.10)
The columns of R are the eigenvectors of P -1 , while the diagonal of P 0 can be computed from the eigenvalues of P -1 . For n = 3, R is the DCM defined by the half axis unit vectors of the ellipsoid forming a right-handed Cartesian frame.

An alternative representation employed for ellipsoid intersections with a strip in section 3.5 is given by describing the ellipsoid surface as a unit ball w under a linear transform

V E(x c , V) = {x|x = x c + Vw, ||w|| ≤ 1} (3.3.11)
See [START_REF] Pope | Algorithms for ellipsoids[END_REF] for useful transformations between different ellipsoid representations.

Observation models

For both the ESMF and the EKF we consider in a first step GNSS code phase and carrier phase observations. We then integrate UWB ranging observations as a first step to robustify the localization scheme against GNSS outages and investigate how this additional information affects filter performance. In the following GNSS and inter-UAS ranging observation models and their respective linearizations required for the ESMF and the EKF are given. Observation error intervals are derived including parametric uncertainty for set membership estimation. All operations are performed in a local NED frame and satellite positions are assumed to have been transformed into this frame.

GNSS differential code phase observations

GNSS receivers are most popularly used for standalone absolute positioning. However, internally they measure the phase shift between an internal copy of the signals emitted by satellites in space and the received signals. This principle is applied to the sinusoidal carrier phase signal as well the pseudo-random number (PRN) binary code modulated on top of it (phase shifts are usually converted to units of distance for convenience).

Due to the principles employed to measure both phase shifts, the measurement noise of code phase observations exceeds the millimeter-level carrier phase noise by several levels of magnitude. Carrier phase observations, on the other hand, are inherently ambiguous due to their periodicity and short wavelength.

From both phase shifts, range information can be derived, making assumptions about the speed of signal propagation through the atmosphere. Variations of the speed of propagation are at once the major source of error (see e.g. [START_REF] Kaplan | Understanding GPS: Principles and Applications[END_REF] for a broader treatment of GNSS ranging principles) and highly correlated for closeby receivers. This spatial correlation is exploited by differential GNSS positioning techniques. By subtracting observations of two receivers, exogenous errors can be almost completely eliminated.

By differencing, absolute position information is lost, and instead the relative position between receivers can be computed. In the following, we employ both differential code phase and carrier phase observations, each to benefit from its respective complementary qualities. GNSS code phase observations, often denoted as pseudo ranges, of UAS i w.r.t. satellite p in units of distance are given by ρ p,i =r p,i + c(∆T i -∆T p ) + ∆r I,p,i + ∆r T,p,i ∆r E,p,i + ∆r M,p,i + ∆r DLL,p,i

with r p,i = ||p i -p p || the slant range between receiver i and satellite p, c the speed of light, ∆T i the receiver clock bias, ∆T p the satellite clock bias, ∆r I,p,i the ionospheric range error, ∆r T,p,i the tropospheric range error, ∆r E,p,i the ephemeris error, ∆r M,p,i the range error due to multipath and ∆r DLL,p,i the code phase tracking tracking noise. Most of these code phase observation errors are highly correlated spatially and can largely be eliminated by forming single differences at epoch k between two receivers on board UAS i and UAS j, leading to

∇ρ p,i,j =ρ p,i -ρ p,j (3.4.2)
=(r p,i -r p,j ) + c(∆T i -∆T j )
+ (∆r M,p,i -∆r M,p,j+ + ∆r DLL,p,i -∆r DLL,p,j )

=∇r p,i,j + c∇∆T i,j + ∇ p,i,j (3.4.3)
where ∇ p,i,j lumps together residual differential observation errors due to uncorrelated hardware noise and multipath. Note in particular that since UAS flying in formation are typically close in altitude, tropospheric propagation delays are highly correlated between receivers and thus largely cancel out. It is convenient to then form double differences w.r.t. to a reference satellite to remove the receiver clock bias common to all observations, thereby later reducing the localization filter state dimension by one. Double differencing single differences w.r.t. a reference satellite r (usually the maximum elevation satellite in view to minimize noise and multipath contamination of double difference observations), leads to

∇ r ρ p,i,j = ||p i -p s || -||p i + d i,j -p s || -||p i -p r || + ||p i + d i,j -p r || + ∇ r p,i,j (3.4.5)
with the relative position vector d i,j = p j -p i between receivers.

Parametric uncertainty

In reality, the first receiver position p i is only known with meter-level standalone precision. Furthermore, broadcast ephemeris are only local approximations of the true satellite orbits, and RMS deviations between predicted and true satellite positions up to tens of meters are recorded [START_REF] Kim | A Long-term Analysis of the GPS Broadcast Orbit and Clock Error Variations[END_REF][START_REF] Warren | Broadcast vs. precise GPS ephemerides: a historical perspective[END_REF], adding additional parametric uncertainty to the observation equation 3.4.5. This parametric uncertainty is traditionally neglected in the GNSS literature [START_REF] Kaplan | Understanding GPS: Principles and Applications[END_REF] due to their submillimeter impact owing to high GNSS satellite orbits (roughly 2 • 10 4 km), i.e its effect tends to "disappear" in the always present receiver hardware observation noise. The first order Taylor expansion of the observation equation (3.4.5) w.r.t. to d about d = 0 leads then to the familiar linear double difference pseudo range observation equation

∇ r ρ p,i,j = H p,r d i,j + R ρ + ∇ r p,i,j (3.4.6)
with the geometry matrix H p,r = ∂∇ r ρ ∂d d=0 The small remainder R ρ is again commonly neglected in the differential GNSS literature. When applying probabilistic filters such as the EKF, neglecting the ensemble of these small effects can be a reasonable simplification due to their millimeter-level magnitude compared to the typically meter-level code phase hardware noise.

Interval inclusion of linearization errors

The deterministic nature of the ESMF requires however guaranteed bounds on the observation error, forcing us to actually quantify these small error terms. To obtain these bounds, we follow in principle the approach proposed in [START_REF] Scholte | A nonlinear set-membership filter for on-line applications[END_REF] and perform an interval evaluation of the Taylor expansion's remainder term. In the following we consider without loss of generality one element y i of an output vector y. The Lagrange remainder being bounded by the interval inclusion

[R i ([x])] = ([x] -x) T ∂ 2 h i ([x]) ∂x 2 ([x] -x) (3.4.7)
interval bounds are obtained by evaluating R i over [x], the interval containing the state. In [START_REF] Scholte | A nonlinear set-membership filter for on-line applications[END_REF], no parametric uncertainty is considered, i.e. it is assumed that the output

y i = h i (x) (3.4.8)
depends only on the state. To deal with observation equations subject to parametric uncertainty of type

y i = h i (x, θ) (3.4.9)
where interval vector bounds on the parameters θ are known, we proceed as follows.

By adopting an augmented state

x = x T θ T (3.4.10)
we can write

y = h(x) (3.4.11)
The state as well as the parameters are bounded by an interval vector [x] with the interval center x0 corresponding to the current state point estimate and nominal parameter values. The first order interval Taylor expansion about x0 is then given by

[y i ] = h(x 0 ) + ∂h ∂ x x=x 0 ([x] -x0 ) + [R i ([x])] (3.4.12)
We split the gradient

[y i ] = h(x 0 ) + ∂h ∂x x=x 0 (x -x 0 ) + ∂h ∂θ θ=θ 0 ([θ] -θ 0 ) + [R i ([x])] (3.4.13)
Evaluating the last two terms of equation (3.4.13) over [θ] and [x] respectively provides then interval bounds on the error of the combined parametric and linearization error. This procedure to incorporate parametric uncertainty is a general one, but depending on the structure of the observation equation, may not provide very tight interval enclosures as we will see in the following. In some cases the degree of overapproximation is acceptable. in other cases, as will be the case for inter UAS range observations, interval evaluations of the parametric error tailored to the problem enable sharper enclosures.

Application to pseudo range double differences

Before applying this method to the double difference observation equation (3.4.5), it is convenient to write it as

∇ρ = 4 i=1 ||E i x|| (3.4.14) = 4 i=1 (x T E T i E i x) 1/2 (3.4.15) = 4 i=1 (x T M i x) 1/2 (3.4.16) with x = (d T p T p T s p T r ) T (3.4.17) x = (d T θ T ) T (3.4.18)
and where the E i matrices form the vector argument of each vector norm term. For instance, for the first term of 3. In this compact form the first and second derivative can be formed conveniently as

∂∇ρ ∂ x = 4 i=1 (x T M i x) -1/2 xT M i (3.4.21) ∂ 2 ∇ρ ∂ x2 = 4 i=1 (x T M i x) -1/2 M T i -M T i x(x T M i x) -3/2 xT M i (3.4.22)
Interval bounds on the observation model error due to linearization and parametric uncertainty can thus be obtained by evaluating

[R ∇ρ ] = ∂∇ρ([x]) ∂θ θ=θ 0 ([θ] -θ 0 ) + ([x] -x0 ) ∂ 2 ∇ρ([x]) ∂ x2 ([x] -x0 ) T (3.4.23)
A sufficiently small degree of overapproximation can thus be obtained, see figure 3.2 for an example. A large number (> 10 3 ) of samples is drawn from an uniform distribution over x to approximate the true shape of the set containing the error terms. Note that the obtained intervals are quite expectedly not sharp due to dependency. The order of magnitude and center of the error term is however well (outer) approximated. Note furthermore that due to the small magnitude of the error term, the achieved tightness of the interval is fully sufficient. Regarding implementation, intervals bounding the receiver standalone position can be found by again applying set membership filtering to the standalone positioning problem, using algorithms such as proposed by [START_REF] Drevelle | High integrity GNSS location zone characterization using interval analysis[END_REF] or by again applying an ESMF. Note that there is a correlation between satellite ephemeris error bounds and receiver standalone position bounds, since ephemeris errors directly affect receiver positioning errors.

No published assessment of maximum broadcast ephemeris errors is known to the author. This constitutes an open problem for set membership filtering. As a workaround, we assume a pessimistic satellite position error interval of 100m in each coordinate.

Time differenced differential carrier phase observations

When basing localization only on pseudorange observations, assumptions need to be made about the maximum displacement of a UAS between samples for the filter propagation step. This in turn requires an accurate model permitting to formulate upper bounds on motion variables such as maximum accelerations, minimum turn radii etc.

In the following we describe how time differenced carrier phase observations can be integrated in a relative position estimation filter that provides accurate tracking by treating position increments as virtual system inputs and decouples the filter propagation step from aircraft motion. Note however that this argument is only valid at the discrete sampling instants, i.e. whenever carrier phase observations are available. In between with r p,i the slant range between receiver i and satellite p c the speed of light, ∆T i the receiver clock bias, ∆T p the satellite clock bias, ∆r I,p,i the ionospheric range error, ∆r T,p,i the tropospheric range error, N p,i the integer ambiguities in carrier signal cycles, ∆r E,p,i the ephemeris error, ∆r M,p,i the range error due to multipath, ∆r P LL,p,i the PLL carrier tracking noise and λ p the carrier signal frequency. Note that the effect of ionospheric delay is of inverse sign compared to pseudo ranges. Just as is the case for code phase observations, forming the difference between observations of two closeby receivers eliminates common range errors and forming double differences removes the differential receiver clock error. Double difference carrier phase observations with the reference satellite r and the UAS i and j are given by ∇ r Φ p,i,j =Φ p,i -Φ p,j -Φ r,i + Φ r,j (3.4.25) =r p,i -r p,j -r r,i + r r,j + λ(N p,i -N p,j -N r,i + N r,j ) (3.4.26)

+ ∇ r Φ p,i,j
where ∇ r Φ p,i,j lumps together residual differential observation errors due to uncorrelated hardware noise and multipath. Using the linearized double difference range observation equation (3.4.5), we can write carrier phase observations as

∇ r Φ p,i,j =H p,r d i,j + R Φ + λ(N p,i -N p,j -N r,i + N r,j ) + ∇ r Φ p,i,j (3.4.27) =H p,r d i,j + R Φ + λ∇ r N p,i,j + ∇ r Φ p,i,j (3.4.28)
where R Φ combines the linearization remainder term and parametric observation model uncertainty due to uncertainty in satellite and receiver position. With 3.4.28 the vector of double difference observations to n satellites is

      ∇ r Φ 1,i,j ∇ r Φ 2,i,j . . . ∇ r Φ n,i,j       =       H 1,r H 2,r . . . H n,r       d i,j +       R Φ,1 R Φ,2 . . . R Φ,n       +       λ∇ r N 1,i,j λ∇ r N 2,i,j . . . λ∇ r N 2,i,j       +       ∇ r Φ 1,i,j ∇ r Φ 2,i,j . . . ∇ r Φ n,i,j       (3.4.29)
∇Φ =Hd + R Φ + ∇N + ∇ Φ (3.4.30)
Note that from here on we drop unambiguous indices for the sake of readability. Computing d i,j from equation (3.4.30) requires first fixing the unknown integer ambiguities e.g. by applying some standard method such as LAMBDA [START_REF] Teunissen | A new method for fast carrier phase ambiguity estimation[END_REF]. Time-differencing observations locally provides another approach that has originally been introduced for precise offline short term trajectory reconstruction [START_REF] Traugott | A timerelative approach for precise positioning with a miniaturized L1 GPS logger[END_REF]. Forming time differences between observations of one receiver removes, in the absence of cycle slips, carrier phase cycle ambiguities or their respective single or double differences. Furthermore, atmospheric delays largely cancel out due to their slow variation over time and resulting high correlation over consecutive epochs. Variations of atmospheric delays over multiple epochs however unavoidably introduce drift in the standalone time-differenced carrier phase solution. In double differences, these error sources are already canceled. We form thus time differenced double difference ( The interval least-squares solution of 3.4.35 is then given by

[∆ d] =(H T k H k ) -1 H T k [∆∇Φ] (3.4.36)
Note that this is a modification of the accumulation strategy studied in [START_REF] Traugott | Precise flight trajectory reconstruction based on time-differential GNSS carrier phase processing[END_REF] and inherits the same attractive properties. Most importantly, the property that only a set of at least four common observations is required over any two consecutive epochs. As a consequence, the availability of the solution 3.4.36 is higher than when using a fixed reference epoch, especially in the dynamic environment of a small UAS. Flight test results with a GPS L1 low-cost receiver [START_REF] Bolting | A GPS Local -Area Augmentation System for small UAS using consumer -grade single frequency receivers[END_REF] suggest that although loss of lock occurs frequently, e.g. in sharp turns, a set of > 4 observations can always be maintained 1 . Note that when using a fixed reference epoch, losing lock to a sufficient majority of the satellites in view between the current and the reference epoch even only once, or a change in the set of satellites in view due to orbital motion can actually break availability of the carrier phase position increment altogether. For any pair of two UAS, the relative position vector interval can thus be propagated in time by

[d k ] = [d k-1 ] + [∆ dk ] (3.4.37)
or alternatively in center-range notation (see 3.3)

[d k ] = [d k-1 ] + ∆ dk + [w ∆d k ] (3.4.38)

Inter-UAS range observations

The range between the antenna centers of UWB ranging devices on board UAS i and j is given by

r i,j,k = ||d i,j,k + R eb,j r b j -R eb,i r b i || + w r,k (3.4.39)
with the positions r b i , r b j of the UWB antenna centers w.r.t. the GNSS antenna centers in the respective body frame and the device hardware measurement noise w r ∈ [w r , w r ]. Equation (3.4.39) is subject to parametric uncertainty due to uncertainty in the attitudes R eb,i , R eb,j of both UAS. In the following we assume that r b i , r b j are small, i.e. the ranging antennas are close to the GNSS antennas, and that interval bounds of their additive effect on the range observation have been determined. We can then lump their effect on the range observation into a combined range error ŵr,k :

r i,j,k = ||d i,j,k || + ŵr,k (3.
[r] =(d T d) 1/2 + ∂r ∂d d=d ([d] -d) + ([d] -d) ∂ 2 r([d]) ∂d 2 ([d] -d) T + [ ŵ] (3.4.42)
with the gradient and Hessian

∂r ∂d = (d T d) -1/2 d T (3.4.43) ∂ 2 r ∂d 2 = (d T d) -1/2 I 3 -d(d T d) -3/2 d T (3.4.44)
Geometrically, the linearized range observation is the distance from a plane tangent to the sphere centered at UAS i at the point of linearization. While a good approximation at large distances, the increasing curvature of the sphere introduces larger linearization errors as UAS get closer to one another, see figure 3.4.

Evaluating the remainder term of (3.4.42) and forming the interval sum with [ ŵr ] provides an interval bound on the total error of the linearized observation equation. Implementing this generic bounding scheme as introduced in equation (3.4.42) leads however to some overapproximation due to dependency in the involved expressions. The larger the linearization error, the less likely it is that the linearized observations can contribute to contracting the guaranteed state set. That being said, it is desirable to find tight bounds on the linearization error to exploit the maximum amount of information captured by range observations. Explicitly evaluating the error between the nonlinear observation equation (3.4.40) and its linearization, exploiting monotonicity, can however provide tight bounds. The explicit linearization error interval is given by T Using the linearized observation equations and their interval error bounds developed in the previous section, we are now ready to apply linear set membership filtering techniques.

[∆r] =(d T d) 1/2 + g[d] -gd -([d] T [d]) 1/2 (3.

Modified Extended Set Membership Filter

The various existing set membership filters employing ellipsoidal set representations (see e.g. [START_REF] Scholte | A nonlinear set-membership filter for on-line applications[END_REF][START_REF] Durieu | Set-membership estimation with the trace criterion made simpler than with the determinant criterion[END_REF]) follow the common recursive structure introduced in section 3.3, each one having a convenient intuitive geometrical interpretation. In a first step, the guaranteed state set Ω(x k,k , Σ k,k ) of the preceding filter step k is propagated using the dynamics model and an ellipsoidal bound Ω(0, Q k ) on the modeling error. This amounts to forming the Minkowski sum of Ω(x k,k , Σ k,k ) and Ω(0, Q k ). This sum generally not being an ellipsoid, an outer ellipsoidal approximation Ω(x k+1,k , Σ k+1,k ) is formed minimizing some measure of size. Algorithms primarily differ in the selection of this measure. The volume, or equivalently the determinant (see [START_REF] Maksarov | Computationally efficient algorithms for state estimation with ellipsoidal approximations[END_REF] for a survey of minimum volume algorithms) as well as the trace of Σ k+,k have been considered. The determinant criterion is reported to potentially produce very thin, long ellipsoids, while the trace criterion, proportional to the sum of squared half axes, tends to lead to better behaved ellipsoids [START_REF] Durieu | Set-membership estimation with the trace criterion made simpler than with the determinant criterion[END_REF].

In a second step, the intersection of Ω(x k+1,k , Σ k+1,k ) and an ellipsoidal approximation of the consistent state set is approximated by a bounding ellipsoid Ω(x k+1,k+1 , Σ k+1,k+1 ), again minimizing the trace or determinant of Σ k+1,k+1 . These two steps are particularly efficient to implement for linear systems, see eg. [START_REF] Durieu | Set-membership estimation with the trace criterion made simpler than with the determinant criterion[END_REF], exploiting firstly the fact that ellipsoids stay ellipsoids under linear transformations and secondly that the set inversion involved in computing the consistent state set can be performed without approximation by matrix inversion.

In the spirit of the Extended Kalman Filter, the Extended Set Membership Filter proposed in [START_REF] Scholte | A nonlinear set-membership filter for on-line applications[END_REF] extends ellipsoidal set membership filtering to nonlinear systems by local linearization. Ellipsoidal bounds on the resulting linearization errors are employed and lumped together with the propagation error bounds and measurement error bounds of the nonlinear model. The algorithm is rather generic as the propagation and update step provide sets of approximating ellipsoids, each parametrized by a scalar parameter. In [START_REF] Scholte | Online nonlinear guaranteed estimation with application to a high performance aircraft[END_REF][START_REF] Scholte | A nonlinear set-membership filter for on-line applications[END_REF] the authors, while mentioning the existence of closed form minimum trace solutions, propose applying numerical optimization to determine the parameter corresponding to the minimum size ellipsoid in each set. Note that the utilization of numerical optimization does not preclude real-time application, since every member of the given sets is guaranteed to contain the propagated set or intersection respectively, so optimization can be stopped anytime. Closed-form solutions are however preferable due to improved efficiency (especially taking into account the still somewhat limited embedded computing resources of UAS) and guaranteed minimization of the selected size criterion. As can be shown, the minimum trace solution presented in [START_REF] Durieu | Set-membership estimation with the trace criterion made simpler than with the determinant criterion[END_REF] can directly be applied to the propagation step. On the other hand, no minimum trace or minimum volume closed form exist solutions for the update step in its form given by [START_REF] Scholte | A nonlinear set-membership filter for on-line applications[END_REF]. In the following we present a modified, fixed time version of the ESMF, adopting the closed form minimum trace ellipsoid approximations proposed in [START_REF] Durieu | Set-membership estimation with the trace criterion made simpler than with the determinant criterion[END_REF]. A further modification we make to the ESMF is the addition of an input vector acting on both the state transition and outputs. This will prove useful for propagating the relative position by time differenced carrier phase observations and to take the position offset of UWB ranging antennas into account.

Algorithm

We provide first the algorithm in general terms, followed by specifics for relative UAS localization.

Propagation The state transition is given by

x k+1,k = f (x k,k ) + u k (3.5.1) Σ k+1,k = A k Σ k+1,k 1 -β k A T k + Q k β k (3.5.2)
where

A = ∂f (x k ) ∂x x=x k,k
. The parameter β k is computed following [START_REF] Durieu | Set-membership estimation with the trace criterion made simpler than with the determinant criterion[END_REF] to minimize the trace of Σ k+1,k to

β k =    tr Σ 1/2 k,k (Σ 1/2 k,k ) T 1/2 + tr Q 1/2 k,k (Q 1/2 k,k ) T 1/2 tr Q 1/2 k,k (Q 1/2 k,k ) T 1/2    -1 (3.5.3)
Update The update of the standard ESMF algorithm is a one-step matrix operation. Although closed-form solutions are mentioned in [START_REF] Scholte | A nonlinear set-membership filter for on-line applications[END_REF] as an option, numerical optimization has to be employed to minimize the trace of the resulting ellipsoid when implementing the one-step update, see [START_REF] Durieu | Set-membership estimation with the trace criterion made simpler than with the determinant criterion[END_REF]. The alternative update step presented in [START_REF] Durieu | Set-membership estimation with the trace criterion made simpler than with the determinant criterion[END_REF] exploits the fact that a closed form minimum trace solution is possible for the intersection of ellipsoids and strips. Each element of the observation vector with its corresponding error interval defines a pair of hyperplanes -a strip -in state space, in turn defining the consistent space set as all states enclosed by the hyperplanes. By consecutive intersection of the propagated state ellipsoid with each strip, an approximation of the minimum trace ellipsoid intersection is obtained. From the linearized observation equation The strips are then recursively intersected with the propagated state ellipsoid Ω(x k+1,k , Σ k+1,k ).

      y k+1,1 y k+1,2 . . . y k+1,m       =       C 1 C 2 . . . C m       x +       w k+1,1 w k+1,2 . . . w k+1,m       (

Recursive intersection

Define

M 0 =Σ -1 k+1,k (3.5.7) c 0 =x k+1,k (3.5.8
)

For i = 1...m
For each element of the observation vector the two hyperplanes are given by Plane 1 :

C i ||C i || -1 x =(y k+1,i -W k+1,i )||C i || -1 =p + Plane 2 : C i ||C i || -1 x =(y k+1,i + W k+1,i )||C i || -1 =p - (3.5.9)
with their common unit normal vector C i ||C i || -1 and the respective distances to the origin p + , p -. As required by the intersection algorithm, we first check for each plane if it intersects the ellipsoid Ω(x k+1,k , Σ k+1,k ), and if not, compute a parallel hyperplane that is tangent to the ellipsoid, see3 appendix A. The resulting tightened hyperplanes being given by Tightened plane 1 : The corresponding tightened strip has the same intersection with Ω(x k+1,k , Σ k+1,k ), but simplifies the intersection algorithm, see [START_REF] Durieu | Set-membership estimation with the trace criterion made simpler than with the determinant criterion[END_REF]. After normalizing with the new error bounds

C i ||C i || -1 x =p + t Tightened plane 2 : C i ||C i || -1 x =p - t ( 3
y k+1,i =W -1 k+1,i ŷk+1,i (3.5.13)
C i =W -1 k+1,i C i (3.5.14)
we obtain a strip

S(y k+1,i , C i ) = {x : |y k+1,i -C i x| ≤ 1} (3.5.15)
in the form required for the intersection algorithm. The intersection of the ith strip and the intermediate ellipsoid

Ω(c i-1 , M -1 i-1
) is then computed as follows, adapted from [START_REF] Durieu | Set-membership estimation with the trace criterion made simpler than with the determinant criterion[END_REF]:

P =M -1 i-1 δ =y k+1,i -C i c i-1 g =C i PC T i γ =C i P 2 C
T i µ = tr(P)

β 1 = 3 g β 2 = g(µ(1 -δ 2 ) -γ) + 2(gµ -γ(1 -δ 2 )) g 2 (gµ -γ) β 3 = µ(1 -δ 2 ) + γ g 2 (gµ -γ) if β 3 > 0 the intersection does not modify Ω(c i-1 , M -1 i-1
) , i.e. both tightened hyperplanes are tangent and

M i = M i-1 c i = c i-1
Otherwise the minimum trace ellipsoid bounding the intersection is computed as

M = 3β 2 -β 2 1 9 N = 9β 1 β 2 -27β 3 -2β 3 1 54 θ = acos N (-M 3 ) 1/2 q =2(-M ) 1/2 cos θ 3 - β 1 3 α * = 1 1 + q M α * =α * M i-1 + (1 -α * )C i C T i c i =M -1 α * α * M i-1 c i-1 + (1 -α * )C i y k+1,i δ α * =c T i-1 M i-1 c i-1 + (1 -α * )y 2 k+1,i -c T i M α * c i M i =(1 -δ α * ) -1 M α * After the last intersection xk+1,k+1 =c m Σ k+1,k+1 =M -1 m
provide the updated guaranteed state ellipsoid at sample k + 1.

Application to relative localization

We now apply the modified ESMF to the problem of cooperative relative localization between UAS. In a most basic configuration, each UAS runs a bank of filters, one for each other UAS. Computational effort could further be reduced by exchanging set estimates between UAS. In the following we consider the case of one UAS tracking the relative position of another one.

To initialize the guaranteed state ellipsoid, the maximum error interval of the differential pseudo range position dρ,0 is computed by interval arithmetic from the least-squares solution of 3.4.6

[ dρ ] =(H T 0 H 0 ) -1 H T 0 [∇ r ρ 0 ] (3.5.16)
The Minimum Volume Ellipsoid (MVE) enclosing this interval then provides the initial guaranteed state set Ω(x 0,0 , Σ 0,0 ).

The state comprises the relative position vector between UAS i and UAS j

x k =d i,j,k (3.5.17)

f (x k ) =d i,j,k + ∆ dk + w ∆d,k+1 (3.

5.18)

A =I 3 (3.5.19) and observations are double differences of pseudo ranges

y k =∇ r ρ k (3.5.20) =H p,r d i,j,k + w ∇ r ρ,k (3.5.21) C =H p,r (3.5.22)
Note that we use the linearized observation equation here due to the small linearization error. The observation error is bounded by the interval sum of the linearization error interval given by equation (3.4.23) and the interval of double difference pseudo range observation noise obtained from zero-baseline observations. Note that in contrast to the standard ESMF, the interval vector bounds on pseudo range observations (see section 3.4.1) can directly by used and do not need to be approximated by a bounding ellipsoid. The propagation error ellipsoid Ω(0, Q k ) is obtained as the MVE of the interval of carrier phase position propagation errors w ∆d,k , see equation (3.4.38).

Incorporating UWB ranging

Range observations are integrated into the filter by augmenting the observation vector to It facilitates an informed judgment of the merits of the proposed set membership filter to compare it to its stochastic EKF counterpart, to which we will consequently turn our attention in the following section.

y k = ∇ r ρ k r i,j,k (3.5.23) = H p,r d k (i, j) + w ∇ r ρ,k ||d i,j,k || + wr,k (3.

Extended Kalman Filter

The Extended Kalman Filter appears to be, along with the Unscented Kalman Filter, the most popular state estimator providing a confidence estimate for nonlinear systems. For the sake of comparison, we implemented two EKFs that are the analogue of the ESMF schemes presented in the previous section. The state transition and observation equations for the EKF are identical to those used for the ESMF, except that stochastic instead of deterministic bounds on the prediction errors and observation errors are assumed. Accordingly, Q KF and R KF are the corresponding process noise and measurement noise covariance matrices.

For the GNSS-only filter we obtain an estimate of the elements of R KF by fitting a normal distribution to zero-baseline double difference samples collected with a low-cost GPS receiver [START_REF] Bolting | Differential GPS for small UAS using consumer-grade single-frequency receivers[END_REF], see figure 3.5. An approximation of R KF is then

R GN SS = σ 2 ∇ r ρ I m (3.6.1)
where σ ∇ r ρ is the standard deviation of the fitted normal distribution and m is the number of pseudo range double differences. When incorporating UWB ranging, we approximate the unknown error distribution of the linearized range observation by setting

σ r = 3 * wr,k 2 , leading to R GN SSU W B = R GN SS 0 0 σ 2 r (3.6.2)
To approximate Q KF , we first compute an estimate R Φ of the covariance matrix of the carrier-phase TDSD observation noise ∆∇ Φ i,j,k again from zero-baseline observations (see figure 3.6). The covariance of w k can then be computed as

(H T i,j,k H i,j,k ) -1 H T i,j,k R Φ (H T i,j,k H i,j,k ) -1 H T i,j,k T (3.6.3)
Note that the resulting covariance matrix is time-varying due to the time-varying geometry matrix H i,j .

Equipped with these covariance matrices, we follow the standard prediction-correction scheme of the EKF.

Derivation of trajectory planning constraints

The ellipsoid produced by the modified ESMF constrains the position vector between GNSS antennas. To obtain an ellipsoid that contains all parts of the UAS, we cover the airframe by another ellipsoid Ω(0, P ac ) centered at the GNSS antenna center and approximate their Minkowski sum by its minimum trace bounding ellipsoid with the same method used in the ESMF propagation step, see equations (3.5.2)-(3.5.3). The result is an ellipsoid

Ω(x, Σ) = Ω(x, Σ k,k 1 -β + P ac β ) (3.7.1)
guaranteed to completely contain the tracked UAS. The ellipsoid Ω(x, Σ) can directly be integrated as a nonlinear constraint into trajectory planning algorithms. It is however convenient to derive a set of constraints linear in the UAS position to enable use of fast LMPC based trajectory planning (see section 2.5.1).

We propose computing an approximate convex hull where the number of faces is a design parameter, representing a tradeoff between tightness of fit and complexity of the resulting polyhedron.

Finding an approximate convex hull

We adopt the algorithm proposed in [START_REF] Xu | An approximate algorithm for computing multidimensional convex hulls[END_REF] for the computation of an approximate convex hull of a set of points in R n . As shown below, it can be applied to find an outer polyhedral approximation of ellipsoids as well. The algorithm is briefly outlined here, see [START_REF] Xu | An approximate algorithm for computing multidimensional convex hulls[END_REF] for a detailed derivation. We first define a set of direction unit vectors n 1 , n 2 ...n N that are evenly spread on the unit sphere. One way to achieve this even sampling is by adopting spherical coordinates

n s i = r Θ φ T (3.7.2)
where r is distance from the origin, Θ the zenith angle and φ the azimuth angle. Given a vector of angles Θ = θ 1 , θ 2 ...θ n evenly spaced on the interval [-π 2 , π 2 ] and φ = φ 1 , φ 2 ...φ n evenly spaced on [0, 2π], a set of unit direction vectors results to

(1 Θ 1 φ 1 ) T , (1 Θ 2 φ 2 ) T , ...(1 Θ n φ n ) T
Memory consumption can be reduced by sampling one half space of the unit sphere and considering each unit vector and its negative. This is left as an implementation detail, as it messes up notation. After converting back to Cartesian coordinates

n i =    sin Θ i cos φ i sin Θ i sin φ i cos Θ i    (3.7.3)
each of these direction vectors n i is the normal vector of a plane given by

n T i d = p i (3.7.4)
where d is a point on the plane.

Application to Ellipsoids

Given the guaranteed state ellipsoid Ω k (x k , Σk ), for each of the direction unit vectors n i we find the p i,t that represents a plane tangent to the ellipsoid and normal to n i by applying the method given in appendix A. Each tangent plane forms one of the planes constituting the approximate convex hull. The resulting set of planes can be written as a polyhedron of the form

Ad ≤ b       n T 1 n T 2 . . . n T n       d ≤       p 1,t p 2,t . . . p n,t       (3.7.5)
Existing algorithms such as [START_REF] Bemporad | Decentralized linear time-varying model predictive control of a formation of unmanned aerial vehicles[END_REF] can then be used to derive convex position constraints from the polyhedron 3.7.5.

Evaluation

To evaluate the proposed localization schemes, we consider two UAS in close proximity. Each UAS is assumed to broadcast its carrier phase and code phase observations as well as to take UWB range measurements. We select a half circle trajectory of relative positions of the follower UAS w.r.t. the predecessor, emulating a follower changing stations inside a formation, a maneuver requiring guaranteed but not overly conservative relative position estimation error bounds for safe but fast execution. Recall that for tight formation flight we want a filter that converges fast to useful, i.e. not overly conservative outer estimates of a set that is guaranteed to contain the relative position.

Figures of Merit

We consider three figures of merit to compare ESMF and EKF as well as to judge the impact of UWB ranging observations. First, the trace of the respective ellipsoid matrices (For the EKF we consider the 3σ level surface ellipsoid), as it gives a scalar measure of the ellipsoid size and thus its relative conservativeness. Second, the distance between the respective ellipsoid center and the relative position vector, representing the usefulness of the ellipsoid center as point state estimate. Third, plugging the true relative position into equation (3.3.9) provides a measure of distance to the ellipsoid surface indicating whether the guaranteed state set indeed contains the true state at all times.

GNSS observation errors

When comparing the EKF, based on a normally distributed, i.e. unbounded process and measurement error model, and the ESMF and BPDMF, based on interval noise bounds, care has to be taken when it comes to simulating observation errors, as the assumed distribution inherently favors one or another filter formulation. While a simulated normal distribution obviously unduly favors the EKF, a uniform distribution is even further from realistic. From inspection of real zero baseline code phase double difference samples, it appears that the noise profile could reasonably be approximated by a biased normal distribution. However, even more realism can be introduced by directly using noise samples from zero-baseline experiments in simulation, thereby removing the need for assumptions about the statistical nature of double difference observation errors. This approach has been adopted here to simulate both differenced carrier phase as well as differenced code phase observations.

UWB ranging observation errors

Our field experiments [START_REF] Lagneaux | Integration d'un capteur de distance dans un système de cinématique temps réel appliqué aux drones[END_REF] with a pair of UWB ranging devices (DecaWave EVK1000 [START_REF] Ye | Experimental impulse radio IEEE 802.15. 4a UWB based wireless sensor localization technology: Characterization, reliability and ranging[END_REF]) indicate a range error interval excluding antenna position offsets of w r ∈ [-0.5 m, 0.5 m] in unobstructed outdoor conditions. Larger error bounds reported in [START_REF] Ruiz | Comparing Ubisense, BeSpoon, and DecaWave UWB Location Systems: Indoor Performance Analysis[END_REF] are most likely due to the considered metalrich industrial indoor environment that is not representative of the conditions expected within a formation of UAS made mostly of composite material. For simulation we draw observation errors from a uniform distribution.

Satellite constellation A static constellation of 7 satellites on a sphere corresponding to GPS orbits, constrained by a minimum elevation of 15°is simulated.

GNSS only

We first run the ESMF and EKF filter formulations with GNSS observations only. Figure 3.7 shows how the ESMF guaranteed state ellipsoid fulfills its primary property and at all epochs contains the true relative position. The ellipsoid contracts drastically over the first few epochs to a quasi steady state in all three axes with a largest half axis of about 6 m. In comparison, the EKF 3σ level set ellipsoid is by roughly one order of magnitude more optimistic, freeing up considerably more maneuvering space. The ellipsoid contains however from about epoch 45 on no longer the true relative position, see figure 3.9, illustrating the fundamental limitation of Kalman filtering for this application. Note that the EKF estimates an inconsistent confidence ellipsoid in spite of differential pseudo range noise being well approximated by a zero mean normal distribution, i.e. although operating under rather favorable conditions. While both the EKF mean state estimate and the ESMF ellipsoid center approximate the true relative position with sub-meter accuracy (see figure 3.9), the EKF provides a smoother estimate. Note that fluctuations in filter execution times (an interpreted Matlab implementation running on a 4-core i7 CPU on Windows 7 OS) can be attributed to the non-real-time execution environment. 12. The superior point estimation accuracy of the EKF could be exploited by using it for guidance purposes while relying on the confidence ellipsoid of the ESMF for guaranteed collision avoidance. The benefit of incorporating ranging observations is more clearly illustrated by figure 3.13, displaying the ratio of the ESMF a posteriori half axes and trace without and with ranging observations. A reduction in size as measured by the ellipsoid trace of roughly 30% can be observed compared to the GNSS only case. As mentioned in [START_REF] Durieu | Set-membership estimation with the trace criterion made simpler than with the determinant criterion[END_REF], incorporating observations recursively in the update step of the modified ESMF can increase conservativeness of the resulting a posteriori ellipsoid. On the other hand, overbounding the observation error interval by an ellipsoid is likely to increase the conservativeness of the standard ESMF update step compared to the recursive intersection that directly operates on the observation error interval.

We consider two measures to study these opposing effects quantitatively. Firstly, when running the modified ESMF, we compute the minimum trace a posteriori ellipsoid by numerical optimization in parallel with the closed-form algorithm, see figure 3.14. Interestingly, during the initial contraction both an improvement in a posteriori ellipsoid size and increased conservativeness can be observed. Differences diminish over time and after filter convergence both methods are roughly on par. Of more consequence is the size of the converged a posteriori ellipsoid. Running the ESMF in both configurations and forming the ratio of the trace of the resulting ellipsoid trajectories, we observe steady state ellipsoids roughly twice as large (as measured by the trace criterion) when using numerical optimization. Note however that a factor of two in the trace size criterion translates into a factor of √ 2 considering lengths of half axes. Still, this snapshot result suggests that if observation errors are defined by intervals, the overapproximation by ellipsoids when applying the standard ESMF leads to significantly more conservative set estimates, on top of the performance penalty of numerical optimization. It suggests also that the effect of ellipsoid bounding of observation error intervals outweighs the effect of recursive strip intersection.

Overview

We have seen in this chapter how ellipsoidal sets that are guaranteed to contain the relative positions of other UAS can be efficiently computed with a modified version of the Extended Set Membership Filter, relying on low-cost GNSS receiver and UWB ranging hardware. The performed simulations suggest that the resulting algorithm is suited for simultaneously tracking multiple UAS in real-time due to its moderate computational cost. We have demonstrated how the resulting ellipsoidal sets can be outer bounded by polyhedrons to render them useful for path planning algorithms requiring linear position constraints.

Conclusion and future work

Tight formation flight is an enabling technology for a large number of range enhancing techniques such as upwash exploitation, aerial recharging and aerial docking between UAS. The corresponding major technological challenges lie in reliably and accurately maintaining relative position in the wake of another aircraft, while guaranteeing collision avoidance. In this work we address both of these challenges. We present important building blocks that are expected to be of help in taking tight formation flight from the experimental stage and make it a safe and routine technique. We have investigated the critical sensitivity of previously proposed continuous time sliding mode guidance techniques w.r.t. two often neglected properties of real-world aircraft: discrete sampling and inner control loop dynamics. We show how these implementability issues can be overcome by extending established discrete sliding mode techniques. The proposed laws are entirely defined by hard input constraints, an advantage over existing DMSC laws. We show that both existing DSMC control as well as the proposed PDSMC control can be posed as special cases of model predictive control To complement the proposed novel guidance laws, we present a design procedure based on model predictive control techniques to generate minimum-time constrained sliding surfaces. By tightening input and position constraints, these sliding surfaces allow for taking dynamic limitations of the UAS and bounded control errors into account simultaneously to guarantee closed-loop constraint satisfaction. We present a set membership localization algorithm that enables for the first time to give safety guarantees when it comes to relative positioning between aircraft. The presented set membership filtering scheme provides both guaranteed regions and point estimates of the relative position between UAS. It can be employed in a complementary fashion parallel to possible existing point localization algorithms such as the Extended Kalman Filter. We provide a derivation of linearized observation equations and the corresponding interval inclusions of modeling error bounds. In extension of the standard ESMF algorithm, we take not only linearization errors, but also parametric errors into account that are commonly neglected in differential GNSS positioning. The algorithm can readily be applied to manned aircraft, where safety guarantees are key to making commercial manned tight formation flight a reality in the future. Incorporating UWB ranging observations is a first step towards GNSS-independent positioning. Guaranteed localization is complemented by the predictive sliding surfaces employed in this work, allowing to take bounded position tracking errors into account. As a distinctive property, all the presented algorithms are readily implementable using the avionics hardware of a typical small low-cost UAS.

Beyond the pressing need for experimental verification, the work presented here has opened up numerous directions for future exciting work.

Future research Simulation

• Available atmospheric turbulence models such as the Dryden model currently used by the FFB benchmark do not allow taking into account the spatial distribution of turbulence in a formation. Due to the very small distances between UAS in formation flight, turbulence acting on each UAS can be expected to be highly correlated. As a result, the effect of atmospheric turbulence is likely to be overestimated by current simulations. Efforts should be made towards a multiple-aircraft turbulence model.

Heterogeneous formations

• In this work, formations of identical UAS have been considered. Heterogeneous formations have distinct advantages for certain missions. Examples are search and rescue missions in cluttered environments, where a group of low-flying small UAS can explore an area more efficiently than a single larger UAS at high altitude. Tight formation flight holds potential in this context for two reasons. First, due to their stronger upwash, larger UAS could improve the endurance of following smaller UAS, leading to better overall formation endurance. Secondly, larger UAS could carry batteries and recharge smaller UAS in air, again improving overall group endurance. The stronger wake of larger UAS will make avoiding certain regions of the wake necessary. The constrained sliding surface presented in this work can readily accommodate wake region constraints.

Scalability

• As is illustrated in section 2.6, rigid formations face scalability problems when it comes to maneuvering. This issue should be treated in a global manner in a future research effort. Two lines of attack appear promising. Firstly, for a given formation, new constraints could be introduced into trajectory planning to ensure indi-vidual constraints satisfaction on the level of each formation member. Secondly, dynamically modifying the formation shape during maneuvering could enable even large formations to fly on trajectories designed to be feasible for one UAS. To give one example, by warping a V-formation into a stacked vertical formation, a horizontal turn feasible for one UAS becomes feasible for a formation of arbitrary size.

PDSMC guidance and load factor control

• The PDSMC guidance law has been implemented on the vertical channel of the FFB benchmark. It is the next step to apply PDSMC control to the coupled vertical/lateral channels, followed by experimental evaluation.

• The proposed PDSMC guidance law relies on sufficiently accurate local prediction models of the load factor and bank angle tracking loops. Adaptive control techniques could improve the proposed guidance schemes in two ways. Firstly, indirect adaptive control techniques, i.e. online learning could be employed to learn better models of the inner loop dynamics. Secondly, adaptive control on the inner loop level could be employed to track reference models corresponding to the PDSMC prediction models. The second option appears more promising, since it alleviates at once the need to accurately identify an aerodynamic model currently required for LQ tracking law synthesis.

• Plasma actuation, while still at an early development stage (about TRL 3 as first wind tunnel [START_REF] Nelson | Modification of the flow structure over a UAV wing for roll control[END_REF] and experimental results for UAS actuation are reported ) allows for very fast flow manipulation [START_REF] Patel | Plasma actuators for hingeless aerodynamic control of an unmanned air vehicle[END_REF]. Once this technology has become more mature, it can remove the main obstacle, slow input dynamics, for the application of continuous time higher order sliding mode control, with its distinct advantages outlined in chapter 2, and continous time HOSM should be revisited for the vertical and possibly the lateral channel.

• The PDSMC optimization problem 2.5.52 could be posed employing a prediction model sampled at the sampling time of the inner loops, imposing the same σ contraction constraint on the last state, enlarging the degrees of freedom of the control. It is considered a worthy future study to determine the degree of performance loss due to the limitation to constant control inputs over the contraction horizon.

• In this work estimated load factors are assumed to be available. Accurate 3D load factor estimation with the consumer-grade MEMS sensors found on todays small low-cost UAS is an interesting problem in itself and deserves future study.

Localization

• In evaluating the set membership localization filter presented in chapter 3, no effort has been made to optimize the trajectory of the follower, i.e. by maximizing some observability measure. Including the contraction of the guaranteed state set into the design of the approach trajectory, based on some kind of observability metric (see e.g. [START_REF] Ponda | Trajectory optimization for target localization using small unmanned aerial vehicles[END_REF] considering the Fisher information matrix) contraction time as well as final size of the guaranteed state set could be systematically minimized. Consider also the QP based approach for bearings-only observability optimized approach trajectories [START_REF] Grzymisch | Optimal rendezvous guidance with enhanced bearings-only observability[END_REF] that appears to be a good candidate to be adapted to the set membership localization scheme considered in this work. The projection matrix G fulfills the usual conditions to ensure asymptotic convergence to the reference state in sliding mode. Applying identical adaptive super twisting controllers [START_REF] Castañeda | Adaptive super twisting flight control-observer for a fixed wing uav[END_REF] to each subsystem and simulating a chain of 5 vehicles, where the first vehicle is excited by a periodic disturbance, we observe two phenomena. Firstly, during the state convergence phase, control inputs grow with the vehicle index. Secondly, while the first subsystem maintains at least first order sliding mode (we check only σ = 0 here) from t = 0 on, the following vehicles cannot maintain sliding mode and enter a reaching phase, exhibiting control inputs growing with vehicle index as well.

.1. 4 )Figure 1 . 1 :

 411 Figure 1.1: Example time series of SVM load factor disturbances

Figure 1 . 2 :

 12 Figure 1.2: Cularis glider (image taken from www.multiplex-rc.de)

Figure 1 . 3 :

 13 Figure 1.3: Normalized differential thrust and normalized aileron deflection of follower UAS i, ∆x = 2b, ∆z = 0, V a = 15 m s , variation of lateral separation

1 Figure 1 . 4 :

 114 Figure 1.4: Normalized differential thrust and normalized aileron deflection of follower UAS i, ∆x = 2b, ∆z = 0, V a = 15 m s , variation of vertical separation

Figure 1 . 5 :

 15 Figure 1.5: Benchmark trajectory, curvature over arc length

3 . 1 )

 31 Cruise flightThe cruise flight primitive adds an additional waypoint to the trajectory in function of a desired incremental trajectory heading ∆ψ and a desired segment

Figure 1 . 6 :

 16 Figure 1.6: Benchmark trajectory: 3D view, projection on x e , y e , altitude profile over arc length

Figure 1 . 7 :

 17 Figure 1.7: Isoview of completed benchmark trajectory in UAV3D

Figure 1 . 8 :

 18 Figure 1.8: N=3 formation in UAV3D

Figure 2 . 2 :

 22 Figure 2.2: Predecessor-follower geometry longitudinal-vertical

Figure 2 . 3 :

 23 Figure 2.3: n z tracking example, Cularis UAS model, DLC LQ control

  .1.5)) can be set to one. The FFB standard vehicle model (see chapter 1.1) is considered, and simulations are run over a sweep of two parameters: inner loop time constant τ ∈ [τ , τ ] = [0.01, 0.5]

Figure 2 . 4 :

 24 Figure 2.4: Mean squared error and maximum vertical position error under sampled super twisting control over envelope of controller sampling time and inner loop tracking time constant, SVM cruise flight scenario with periodic predecessor motion. White areas mark performance indices exceeding one wingspan.

Figure 2 . 5 :

 25 Figure 2.5: Tracking error and sliding variable under sampled super twisting control for corner points of sampling time and inner loop time constant envelope, SVM cruise flight scenario with periodic predecessor motion

Figure 2 . 6 :

 26 Figure 2.6: FFB benchmark, state convergence phase with large initial relative position error in vertical channel, turbulent air, nominal state constraints ∆v g z ∈ [-2, 2]

Figure 2 . 7 :

 27 Figure 2.7: FFB benchmark, state convergence phase with large initial relative position error in vertical channel, turbulent air, tightened state constraints

.5. 75 )Figure 2 . 8 :

 7528 Figure 2.8: Disturbance estimation, T s = 10 -2 s

Figure 2 . 9 :

 29 Figure 2.9: Disturbance estimation, T s = 10 -1 s

Figure 2 . 10 :

 210 Figure 2.10: Sensitivity of disturbance estimation error w.r.t observation sampling time

Figure 2 . 11 :Figure 2 . 12 :Figure 2 . 13 :

 211212213 Figure 2.11: Mean squared error and maximum vertical position error under PDSMC control over grid of controller sampling time and inner loop tracking time constant, SVM cruise flight scenario with periodic predecessor motion

Figure 2 . 14 :Figure 2 .Figure 2 .

 21422 Figure 2.14: Baseline control, T s = 10 -1 s, turbulent air, (a) entire trajectory, (b) closeup of cruise flight phase and first maneuver. The WMS error applies to the cruise flight segment between 15 s and 55 s

. 22 .Figure 2 . 17 :Figure 2 . 18 :Figure 2 . 19 :Figure 2 .Figure 2 . 21 :Figure 2 . 22 :

 222172182192221222 Figure 2.17: PDSMC, T s = 10 -1 s, turbulent air, disturbance estimation, large initial deviation, (a) entire trajectory, (b) closeup of cruise flight phase and first maneuver

D 1, 3 D 1 Figure 3 . 1 :

 3131 Figure 3.1: Relative localization problem in 2D: guaranteed relative position sets and possible ellipsoidal and interval outer approximations

4 . 5

 45 ||p -p s || =||E 1 x|| (3.4.19) =|| 0 3 I 3 -I 3 0 3 (d T p T p T s p T r ) T ||(3.4.20)

Figure 3 . 2 :

 32 Figure 3.2: Interval bounds on pseudorange double difference linearization error including parametric error for two satellites, Monte-Carlo sampling approximating true error set

4 . 40 )

 440 with ŵr ∈ [ ŵr , ŵr ] Linearization To linearize equation (3.4.39) for the ESMF and EKF, we first write range observations as r =(d T d) 1/2 (3.4.41) dropping indices and subscripts since we always consider UAS i and j at sampling instant k. The first order interval Taylor expansion of equation (3.4.41) follows in the form familiar from pseudo range observations

4 . 45 ) 47 )Figure 3 . 4 :

 4454734 Figure 3.4: Range observation linearization error interval bounds and Monte Carlo samples for two separation vector intervals: [d a ] = [90, 110] [-10, 10] [-10, 10] T and [d b ] = [0, 20] [-10, 10] [-10, 10]

  3.5.4) where C i are the rows of C = ∂h(x k ) ∂x x=x k+1,k and w k+1,i ∈ [-W k+1,i , W k+1,i ] 2 we obtain m pairs of hyperplanes defining m strips in R n . If known inputs affect outputs in the known form ỹk =h(x) + g(u k ) + w k (3.5.5) we define a new output y k = ỹk -g(u k ) allowing us again to write y k =h(x) + w k (3.5.6) and to linearize without modifications to obtain equation (3.5.4).

Figure 3 . 5 :Figure 3 . 6 :

 3536 Figure 3.5: Time series and fitted normal distribution of GPS C/A pseudo range double difference noise obtained from zero-baseline measurements (ublox LEA-6T receivers)

Figure 3 . 7 :

 37 Figure 3.7: ESMF, GNSS only: half axes of guaranteed state set, ellipsoid membership of true state, guaranteed state set at last epoch with pseudo range positions and ESMF position estimates

Figure 3 . 8 :Figure 3 . 9 :

 3839 Figure 3.8: EKF, GNSS only: half axes of ellipsoid Σ 3σ , ellipsoid membership of true state for Σ 3σ , Σ 3σ at last epoch with pseudo range positions and KF position estimates

Figure 3 . 3 . 8 . 2

 3382 Figure 3.10: ESMF, GNSS+UWB ranging: half axes of guaranteed state set, ellipsoid membership of true state, guaranteed state set at last epoch with pseudo range positions and ESMF position estimates

Figure 3 .Figure 3 . 12 :Figure 3 . 13 : 3 . 8 . 3 FixedFigure 3 . 14 :

 3312313383314 Figure 3.11: EKF, GNSS + UWB ranging: half axes and trace of ellipsoid Σ 3σ , membership of true state in Ω(x k,k , Σ 3σ,k,k , Σ 3σ at last epoch with pseudo range positions (blue dots) and EKF position estimates (red line)

Figure B. 1 :

 1 Figure B.1: Chain of 5 vehicles, identical sliding surfaces, T s = 10 -4 s

  The depicted four cases represent the corner points of the sampling time and inner loop time constant envelope. Note that these simulation results have to be considered with caution quantitatively due to the severely simplified nature of the SVM, but qualitatively demonstrate the rapid deterioration of TSST performance in the presence of moderate sampling times and moderately fast inner loop dynamics.

[START_REF] Stolle | A vision-based flight guidance and navigation system for autonomous cross-country soaring UAVs[END_REF]

, good tracking performance can be achieved for conditions close to the design model. Performance is however very sensitive w.r.t. both parameters and exceeds the defined acceptable interval for a large part of the parameter space. At moderate sampling rates and with realistic inner loop constants, heavy chattering occurs, see time series of normalized position error and the sliding variable in figure 2.5.

  Consider the error dynamics of a pair of SVMs

	∆ ṗg i (t) =∆v g i (t) + ξ v	(2.5.63)
	=∆v g i (t) + Ω g (t)∆p g	(2.5.64)
	∆ vg i (t) =gn g i+1 + ξ a	(2.5.65)

  [START_REF] Weimerskirch | Energy saving in flight formation[END_REF] to the same SVM tracking problem employed to evaluate TSST performance, see figure 2.4. No measurement noise is considered. PDSMC control shows improvements w.r.t. both performance indices for small sampling times and fast inner loops as well as the desired graceful performance degradation for slower inner loops and lower sampling times. See figure2.13 for a closeup of the reaching phase and a short period of the quasisliding phase under super twisting control and PDSM control under discrete sampling. We have selected for this illustrative simulation a setting where super twisting control

	0.1
	0.08	DDO noisy	DDO noise free	SMO noisy	SMO noise free
	0.06
	0.04
	0.02
	0
	0.25
	0.2	DDO noisy	DDO noise free	SMO noisy	SMO noise free
	0.15
	0.1
	0.05
	0
	0.01	0.02	0.03	0.04	0.05	0.06	0.07	0.08	0.09	0.1

starts to exhibit significant, but still acceptable chattering. Under PDSM control, on the other hand, after a moderate overshoot over the sliding surface no visible chattering is encountered, see also figure 2.12. Encouraged by these preliminary results, we apply the disturbance observer based PDSMC guidance law to the vertical channel of the

  The basic PDSMC scheme (see figure 2.19) again displays maximum tracking errors inferior to the baseline scheme (see figure 2.18). Disturbance estimation improves performance even more visibly than under turbulence, see. figure 2.20. Overall, the tracking errors observed, if maintained in flight tests, are sufficiently small to enable considerable upwash gains. The performance gain is more clearly visible by a direct overlay, see figure 2.21.

  This aspect is an important future extension of the work presented here. While the noise level of pseudo range double differences translates into position errors of the order of meters, carrier phase observations have complimentary properties as they feature a very low, i.e. millimeter level noise component (see figures 3.5 and 3.6), but are biased by unknown integer phase ambiguities due to the periodicity of the carrier signal. GNSS carrier phase observations of UAS i w.r.t. satellite p in units of distance are given by Φ p,i =r p,i + c(∆T i -∆T p ) -∆r I,p,i + ∆r T,p,i λN p,i + ∆r E,p,i + ∆r M,p,i + ∆r P LL,p,i (3.4.24)

  T3D) observations at epoch k w.r.t. an epoch k 0 ∆∇Φ =H k ∇d k -H k 0 d k 0 + ∆R Φ,k + ∆∇ Φ (3.4.31) Figure 3.3: Interval range of two matrix elements of H, approximation of true result set By choosing k 0 = k -1, satellite-receiver geometry varies only very little due to the low airspeeds of small UAS, and we can approximate ∆∇Φ ≈H k ∆d + ∆R Φ,k + ∆∇ Φ (3.4.32) and compute the least squares solution of ∆d = d k -d k0 . To find a bound on the approximation error we consider bounds on the variation ∆H of the geometry matrix H between two consecutive epochs by first writing ∆∇Φ =H k d k -(H k + ∆H)(d k + ∆d) + ∆R Φ,k + ∆∇ Φ (3.4.33) =H k ∆d + ∆H∆d k + ∆H k ∆d + ∆R Φ,k + ∆∇ Φ (3.4.34) By interval arithmetic we can then find interval bounds on the second and third term of equation (3.4.34), see figure 3.3 for an example of the obtained intermediate interval of two elements of ∆H. Note that the inclusion is not tight, but predicts well the nonzero center and order of magnitude of the true range, approximated by a large number of point samples. The interval sum of the last four terms of equation (3.4.34) provides then an outer bound [w Φ ] for the overall observation error of linearized time differenced carrier phase double differences and the observation interval is

	10 -6
	6
	4
	2
	0
	-2
	-4	-2	0	2	4
					10 -6
	[∆∇Φ] =H k ∆d + [w Φ ]	(3.4.35)

By interface variables we mean those aircraft states that are tracked by inner controllers and serve as virtual control inputs for the outer controllers

The matching condition formalizes the property that disturbances enter the system by the same channels as control inputs, in principle enabling the latter to instantaneously eliminate the former

What is more, the ongoing incorporation of concurrent GNSS systems (GPS, Glonass, Compass, Galileo) in low-cost receivers can be expected to improve availability drastically in the near future.

Note that every asymmetric observation error interval can be converted into a new, centered observation and a corresponding symmetric error interval

Details on computing these tightened hyperplanes might also be given in[START_REF] Belforte | An improved parameter identification algorithm for signals with unknown but bounded errors[END_REF], where tightening has first been introduced. Since[START_REF] Belforte | An improved parameter identification algorithm for signals with unknown but bounded errors[END_REF] appears to be available only in print, we provide a derivation in the appendix for the sake of accessibility

Appendix A

Computing tangent hyperplanes of ellipsoids

Given an ellipsoid and a plane E =Ω(x c , Σ) n T x =p with the plane normal vector n, we check first whether the plane intersects the ellipsoid and, if not, find a modified p t that represents a plane parallel to the original one and tangent to the ellipsoid. To considerably simplify doing so, we first apply a linear transform to the problem space that distorts the ellipsoid back into a unit sphere, find the tangent plane in the distorted space and transform it back to the original space. The ellipsoid being a unit ball under a linear transform, the inverse transform is given by

The transformed ellipsoid becomes the unit sphere

where () denotes a vector in the transformed space, and the plane corresponding to the unit vector n in transformed space is given by

APPENDIX A. COMPUTING TANGENT HYPERPLANES OF ELLIPSOIDS

In the transformed space the distance of the initial plane from the surface of the sphere is straightforwardly given by

By replacing p by an appropriately modified p t the plane is made tangent to the sphere. The tangent plane is then transformed back to the original space by T -1

Appendix B Disturbance propagation under local SMC

String stability and its 3D generalization mesh stability are a feature of a three-dimensional formation of vehicles that allow separation errors to stay locally contained. It is a well known fact ( [START_REF] Pant | Mesh stability of look-ahead interconnected systems[END_REF]) that formations under linear control with local feedback are mesh unstable.

Very few solutions based on nonlinear techniques such as Port-Hamiltonian systems, [START_REF] Knorn | Passivity-based control for multi-vehicle systems subject to string constraints[END_REF] and sliding mode control, [START_REF] Bolting | Discrete sliding mode control of small UAS in tight formation flight under information constraints[END_REF][START_REF] Galzi | UAV formations control using high order sliding modes[END_REF], have been proposed so far to tackle this problem and no experimental results are known to the authors. Furthermore, the existing approaches have been applied to very simplified models, basically point-mass approximations, neglecting a large part of the complex dynamics and constraints of a fixedwing UAS. While in ideal sliding mode, the position error dynamics are by definition confined to the sliding surface, independently of adjacent separation errors, implying mesh stability if the system can be kept in ideal sliding mode. It is this insensitivity towards bounded disturbances that makes sliding mode control a very interesting candidate for mesh stable formation control with local state information.

In the following we show that even under conditions allowing for ideal sliding mode (continuous time system, no input dynamics), bounds on control inputs required to maintain motion on the sliding surface grow downstream, making SMC unsuitable for large formations with local feedback.

1D mesh stability in ideal sliding mode In ideal sliding mode after the reaching phase (t > t reach ), the system states are confined to the sliding surface and exogenous disturbances completely compensated for, given that they remain within the design bounds. Consider the simplified vertical motion of a pair of two UAS i and i + 1 in the sliding phase (σ = 0), where

then σ = 0 provides the equivalent control applied by the predecessor i + 1

By the triangle inequality we obtain an upper bound for the equivalent control (the acceleration of the follower vehicle, a i+1 required to completely reject disturbances and maintain sliding motion) as

where ||v|| ∞ denotes max t≥0 (||v(t)||). The first term on the right side of equation (B.0.5) represents the input component compensating disturbances while the second term drives the compensated dynamics along the sliding surface. The second component represents the control input generated by the nominal controller, since without disturbances integral sliding mode is already induced by the nominal controller. In ideal sliding mode, the system states can be made to reach the origin and stay confined to it in finite time by the right choice of sliding surface , i.e. ∆p i (t) = ∆v i (t) = 0 for t > t * > t reach . We denote the time interval t * > t > t reach the state convergence phase. Note that, during the state convergence phase, upper bounds on the accelerations the follower has to exert to maintain sliding motion exceed those of the predecessor. Since the acceleration of the follower enters the next system downstream as disturbance, by recursion this implies that norm bounds on the acceleration demands in a chain of vehicles in sliding mode grow with the vehicle index during the state convergence phase. In quasi sliding mode, the system stays in a reaching phase inside the boundary layer at all times, therefore this issue affects quasi sliding mode in an even more severe fashion. We illustrate this issue by a simple simulation example. Consider a chain of double integrator UAS models. An integral sliding surface designed by applying LQ control results to σ = G x -