AHRS Attitude Heading Reference System ALOFT Autonomous Locator of Thermals

Small fixed-wing Unmanned Aerial Vehicles (UAVs) provide utility to research, military, and industrial sectors at comparably reasonable cost, but still suffer from both limited operational ranges and payload capacities. Thermal soaring flight for UAVs offers a significant potential to reduce the energy consumption. However, without remote sensing of updrafts, a glider UAV can only benefit from an updraft when hitting it by chance.

In this thesis, a new framework for autonomous cross-country soaring is elaborated, enabling a glider UAV to visually localize sub-cumulus thermal updrafts and to efficiently gain energy from them.

Relying on the Unscented Kalman Filter (UKF), a monocular vision-based method is established, for remotely estimating sub-cumulus updraft parameters. Its capability of providing convergent and consistent state estimates is assessed by Monte Carlo Simulations. Model uncertainties, image processing noise, and poor observer trajectories can degrade the estimated updraft parameters. Therefore, a second focus of this thesis is the design of a robust probabilistic path planner for map-based autonomous cross-country soaring. The proposed path planner balances between the flight time and the risk of an outlanding by taking into account the estimation uncertainties in the decision making process. The suggested updraft estimation and path planning algorithms are jointly assessed in a 6 Degrees of Freedom (DOF) simulator, highlighting significant performance improvements with respect to state of the art approaches in autonomous cross-country soaring while it is also shown that the path planner is implementable on a low-cost computer platform. i the holidays, the dinners, the conversations. I am grateful for being your brother.

To my parents: I am grateful for your

But those who trust in the LORD will find new strength. They will soar high on wings like eagles. They will run and not grow weary. They will walk and not faint.

-Isaiah 40:31 (NIV)

To My Parents -With Love and Gratitude, and above all, To the almighty God. This thesis was not a solo effort. I am grateful for every individual who supported me.

First I would like to thank Carsten Döll for providing every kind of help from the very early stage of my project. Thanks a lot Carsten for having been supportive with the scholarship and the valuable exchange we had talking about flight control and soaring techniques. I would also like to thank Yoko Watanabe for the countless hours she dedicated to my work.

Thank you for having shared with me your knowledge and for having helped with improving and concretizing many of the ideas of this work. I am also grateful towards Caroline Chanel, having spent so much time with me talking about path planning concepts. I am also grateful for your encouragement.

My committee, Robert Luckner, Isabelle Fantoni, and Yann Le Gorrec, thank you for your interest in my work, and for the enriching comments and the advice. I would also like to thank the Multiplex Modellsport GmbH & Co.KG for having provided fully equipped Cularis RC gliders that were valuable for collecting flight data and video records.

To my friend Jan Bolting, thank you a lot for the past years. I'm still thinking of the times where we started working together in the TU Braunschweig UAV lab, especially the UAV competitions in Florida, and of course our excursions to Hossegor.

To Leandro Lustosa, I am very happy having met you at the countless ISAE lunch sessions and grateful for having shared an office with you during the last period of my work.

To Timothy van Oorschot, I had never such a nice roommate who turned out to be a divine cook, a balanced and grate character, and a good friend. Thank you also for the unforgettable Atlantic Ocean weekends. I must say these Cap Ferret days were part of my best times in France.

To François Bachoc, thank you for having been a supportive roommate and for your countless advices.

To my brother and sister Simon and Charlotte. Thank you for the phone calls, the reunions,

List of Symbols and Acronyms

If not stated otherwise, scalars are noted as non-bold, vectors as bold, and matrices as bold and/or non-bold uppercase.

Symbols

Introduction

Vision-based autonomous soaring is an unexplored and promising way to extend both the endurance and range of fixed-wing UAVs. This subject poses challenges in atmospheric modeling, flight control, image processing, non-linear estimation and path planning. Among these, the underlying work focuses on the estimation and path planning parts and results in two main contributions:

The first one refers to the vision-based remote estimation of sub-cumulus updrafts. It provides an estimator for remotely sensing updraft parameters that are necessary to perform efficient autonomous cross-country soaring flight.

The second contribution is a robust, map-based path planning algorithm for autonomous cross-country soaring. It balances between the risk of outlandings and the minimization of flight time by considering uncertainties of updraft map parameters.

Background and Motivation

In the past decades, UAVs have proven to be of great utility, facilitating a variety of tasks while simultaneously reducing operating costs. Many applications such as border patrolling or pipeline monitoring require the UAV to cross significant distances, or simply to stay aloft for long periods of time, e.g. remote sensing or atmospheric research. Both these requirements however, challenge especially small and electrically powered vehicles. This is primarily due to one important aerodynamic performance metric, the maximum lift-to-drag ratio. High liftto-drag ratios call for high wing aspect ratios, conflicting with the small size of these vehicles.

Additionally, the small carriage capacity is often utilized in favor of the payload, rather than the battery.

While it has become common practice to interrupt the mission in order to change batteries, the development of batteries with higher energy densities might be a long-term solution. Other solutions aim at weight reduction by using advanced materials such as carbon fiber, or include the design of aerodynamically more efficient airframes. Even the development of flight control laws for UAV formation flight is considered to reduce energy consumption [START_REF] Bolting | Discrete sliding mode control of small UAS in tight formation flight under information constraints[END_REF], mimicking the behavior of migratory birds.

In comparison, both range and endurance of a UAV can also be extended by making use of external energy -provided by "mother nature". For example, efforts were conducted to equip UAVs with solar panels for an in-flight battery recharge [START_REF] Morton | Solar powered UAV: Design and experiments[END_REF][START_REF] North | Autonomous solar UAV for sustainable flights[END_REF]. Another option is to use external energy in order to perform unpropelled flight. Unpropelled flight, widely known as soaring, is complementary to the aforementioned approaches. Compared to the development of high performance batteries, the use of solar panels or advanced materials, soaring solely necessitates suited algorithms and in some cases additional sensors. Most generally, soaring flight comprises the two fields of dynamic and static soaring.

Dynamic Soaring

Horizontal wind gradients are used to increase the kinetic energy in dynamic soaring, a method that is commonly adopted by maritime birds such as Albatrosses. It enables them to stay aloft for hours and to travel significant distances without taking a rest. UAVs sounds promising, a successful realization of this method relies on sensors with an accuracy that enables flight in extreme ground proximity. For maritime applications, the sensor has to be accurate enough to allow a reliable "below"-surface operation [START_REF] Deittert | Engineless unmanned aerial vehicle propulsion by dynamic soaring[END_REF].

Figure 1.1: Flight trajectory of an Albatross bird in dynamic soaring [START_REF] Richardson | Upwind dynamic soaring of albatrosses and UAVs[END_REF] Static Soaring Contrary to dynamic soaring, static soaring necessitates the presence of vertical winds. As part of static soaring, slope soaring, also known as ridge soaring, relies on vertical winds that are generated when horizontal winds hit an obstacle. The procedure is illustrated in Figure 1.2. These objects are typically ridges or steep coasts. In urban environments, obstacles can also be buildings. Numerous studies were presented for generating slope soaring trajectories for UAVs (e.g. [START_REF] Langelaan | Long distance/duration trajectory optimization for small UAVs[END_REF][START_REF] Cutler | Energy harvesting and mission effectiveness for small unmanned aircraft[END_REF][START_REF] White | A feasibility study of micro air vehicles soaring tall buildings[END_REF]). Vertical winds can also be generated by solar radiation, heating the ground and creating columns of rising air. These columns are referred to as thermal updrafts. Thermal soaring for UAVs was first investigated in [START_REF] Wharington | Autonomous control of soaring aircraft by reinforcement learning[END_REF]. As reported by the author of [START_REF] Michael | Autonomous soaring for improved endurance of a small uninhabited air vehicle[END_REF], UAVs with a nominal endurance of two hours can increase the flight time by 12 hours in the summer and 8 hours in the winter respectively -only relying on thermal updrafts. This significant potential has motivated researchers to work on autonomous thermal soaring. The first successful realization of an autonomous soaring system was conducted by the NASA in 2007 [2] using a small-scale glider, as illustrated in Figure 1.3. Its successful flight test campaign has since inspired a lot of succeeding research, for example path and motion planning for autonomous gliders [START_REF] Kahveci | Adaptive LQ control with antiwindup augmentation to optimize UAV performance in autonomous soaring applications[END_REF][START_REF] Nguyen | Energy-constrained motion planning for information gathering with autonomous aerial soaring[END_REF],

updraft estimation [START_REF] Stolle | A vision-based flight guidance and navigation system for autonomous cross-country soaring UAVs[END_REF][START_REF] Stolle | A sigma-point Kalman filter for remote sensing of updrafts in autonomous soaring[END_REF], and cooperative soaring [START_REF] Cobano | Thermal detection and generation of collision-free trajectories for cooperative soaring UAVs[END_REF].

Figure 1.3: NASA Autonomous Soaring UAV [2] Compared to other energy sources such as updrafts caused by ridges, thermal updrafts repeatedly emerge over large geographical regions. This fact is exploited in cross-country soaring where gliders/birds fly beyond the gliding distance from the initial take-off point and perform waypoint navigation. The first and sole realization of an autonomous cross-country soaring flight was published by Dan J. Edwards et al., the authors of [3].

Figure 1.4: ALOFT autonomous glider [3] The work led to the participation in a cross-country soaring challenge for remotely piloted gliders. Figure 1.5: Flight log of the ALOFT glider on day one of the competition [3] With no a priori information on thermal locations in the far environment, the UAV flight path was defined as the direct line between two consecutive waypoints. The vehicle's flight control mode was set to thermal centering mode, when encountering strong enough thermals. Thermals were detected by Inertial Measurement Unit (IMU), Global Positioning System (GPS) and airdata measurements. These quantities can be transformed in order to monitor the vehicle's total energy rate. While the 50 km flight was a giant step towards autonomous cross-country soaring, the application of the system is restricted to weather conditions with high updraft densities. Additionally, the lack of future updraft knowledge results in reactive and therefore inefficient flight trajectories. As a consequence, to this day, autonomous cross-country soaring still bares a significant risk of mission failure and therefore finds only limited application.

In contrast, human glider pilots rely on their vision when performing soaring. More specifically, pilots remotely estimate parameters of sub-cumulus updrafts by observing cumulus clouds that lie in the vehicle's environment [START_REF] Reichmann | Cross-Country Soaring[END_REF][START_REF] Pagen | Understanding the sky[END_REF] (see Figure 1.6).

A Framework for Autonomous Cross-Country Soaring

When solar radiation hits the ground, the air in the vicinity of the earth's surface is heated and therefore expands. The expanding, warm, and moist air rises, and forms a column of lifting air. While rising, the air cools down until it reaches the Cumulus Condensation Level (CCL)

where cumulus clouds form. As long as the updraft remains active, its capping cloud grows.

Once the updraft has vanished, the cloud size begins to decay. The cumulus cloud thus not only indicates the position of the updraft, but its geometry even provides information on the activity of the updraft [START_REF] Reichmann | Cross-Country Soaring[END_REF].

In order to stay aloft by benefiting from these energy sources, autonomous soaring systems necessitate several components. Of these, the most fundamental components are both energystate estimation and thermal centering control. Given the rate of change of the vehicle's energy, bank angle control commands can be generated, for centering the glider within the updraft and thereby increasing its potential energy level.

A more complex soaring task is cross-country soaring, as illustrated in Figure 1.8. Beyond the

A Framework for Autonomous Cross-Country Soaring

elementary task of staying aloft, the vehicle is additionally challenged to cross a distance in order to reach a target destination. To fulfill this task, a series of updrafts has to be visited, before the destination lies in the gliding range of the vehicle. For increased reliability, a more Figure 1.8: Cross-country flight in cumulus soaring conditions advanced system is necessitated enabling soaring flight performances that are comparable to those of a self-propelled UAV. Such a system not only relies on energy-state estimation and thermal centering control, but requires two more components. When aiming to prevent both reactive and "luck"-dependent trajectories, remote updraft estimation functionality is essential. Given a map with estimated updrafts, a suited path planner is additionally required, for guiding the vehicle to its destination while minimizing/maximizing given performance criteria such as flight time.

To this date, such an advanced system does not exist. Therefore, a framework as illustrated in Figure 1.9 is proposed in this thesis. This vision-based framework includes the mentioned system components. More specifically, both a Red Green Blue (RGB) camera and an image processing algorithm form a vision unit that provides vision measurements of cumulus clouds.

Vision, position, and attitude measurements, as provided by the Attitude Heading Reference System (AHRS) and the GPS are fused for remotely estimating positions p u , vanish times t v and strengths w u,rz of the updrafts. The generated updraft map serves as input to a path planner that computes a sequence of updrafts-to-visit, necessitated to reach the destination, as given by the mission. High-level control inputs are transmitted to a flight guidance system that generates bank angles and airspeed commands which are tracked by a set of dedicated flight control laws.

Mission

Related Work

The primary components that enable Vision-Based Autonomous Cross-Country Soaring (VBACCS) are updraft estimation, path planning, and thermal centering control. The following section provides an overview of the state of the art on these fields.

Atmospheric Modeling

The design of autonomous soaring algorithms depends on atmospheric models that describe the meteorological phenomena a glider encounters during soaring flight. Evidently, methods like Large Eddy Simulation (LES), commonly used by meteorologists [START_REF] Golaz | A Large-Eddy Simulation study of cumulus clouds over land and sensitivity to soil moisture[END_REF], provide a rather realistic description of the atmosphere. As these models however tend to be computationally costly, they do not provide a practical solution for the design of soaring algorithms. A measurement based updraft model, specifically dedicated to the design of soaring algorithms and widely-used in the community, is provided in [START_REF] Allen | Updraft model for development of autonomous soaring uninhabited air vehicles[END_REF]. A rather simple model that approximates updrafts as a Gaussian lift distribution can be found in [START_REF] Wharington | Autonomous control of soaring aircraft by reinforcement learning[END_REF]. Both these models are limited in the way that they do not consider horizontal winds. Therefore, Lawrance [START_REF] Lawrance | Autonomous soaring flight for unmanned aerial vehicles[END_REF] has made efforts to provide models that include the interaction between horizontal wind and updrafts.

For vision-based autonomous soaring, not only updrafts but also clouds have to be modeled.

Ideally, these models are backed by measurements of local pressure, humidity, and wind-flow, as well as cloud width measurements. Therefore, meteorologists have equipped small UAVs with sensors, in order to collect in-cloud measurements at high resolution [START_REF] Renzaglia | Monitoring the evolution of clouds with UAVs[END_REF]. While in-cloud data collection is feasible, it remains challenging to collect cloud width measurements over time. Therefore, suited image processing algorithms as well as special observer flight trajectories are required, when aiming to rely on UAV systems. Cloud width measurements can also be collected by an observer that is located on the ground. This was done in [4]. To the author's knowledge, this is the only publication of cloud width measurement data. In this work, these measurements are utilized to derive a simplified updraft cloud model which enables the design of vision-based autonomous soaring algorithms.

Updraft Estimation

The updraft estimation problem can be grouped in two fields: In updraft and remote updraft estimation.

In updraft estimation deals with the estimation of the vehicle's vertical velocity. Therefore, the vehicle's energy-state has to be determined. The authors of [2] propose a filter for estimating the vehicle's energy-state and its derivative from static and total pressure measurements while climbing within the updraft. The position, radius, and strength of the thermals are estimated using a history of total energy rate and aircraft positions, stored in a first-in-first-out queue.

In [START_REF] Daugherty | Improving Autonomous Soaring via Energy State Estimation and Extremum Seeking Control[END_REF], a Savitzky-Golay filter is suggested for estimating the total energy, rate of change of total energy, and the second derivative of the total energy of the vehicle by making use of airspeed and altitude measurements. A recursive in-updraft estimation method, based on the Unscented Kalman Filter (UKF) is presented in [START_REF] Hazard | Unscented Kalman filtering for real-time atmospheric thermal tracking[END_REF], aiming to reduce the computational requirements for the updraft position estimation. Yet another method is established by the author of [START_REF] Lawrance | Autonomous exploration of a wind field with a gliding aircraft[END_REF], where a system for simultaneous exploration and exploitation of static wind fields is suggested, relying on model-free Gaussian Process regression.

The energy state estimation method from [2] has successfully been flight tested and demonstrated how endurance can be increased by relying on energy provided by thermal updrafts.

More specifically, the method was flight tested by the NASA in 17 flights during which a total of 23 thermal updrafts were encountered and soaring flights of up to 40 min were performed.

Similar soaring flight durations of about 40 min are reported by the authors of [START_REF] Andersson | Autonomous soaring: Flight test results of a thermal centering controller[END_REF], which rely on the same energy state estimator to generate thermal centering control inputs.

As stated by the authors of [3,[START_REF] Lawrance | Autonomous soaring flight for unmanned aerial vehicles[END_REF], soaring performances can greatly be improved with in-flight remote sensing of updrafts. Researchers have therefore focused on cooperative [START_REF] Acevedo | Persistent monitoring with a team of autonomous gliders using static soaring[END_REF][START_REF] Andersson | Cooperating UAVs using thermal lift to extend endurance[END_REF], and distributed autonomous soaring [START_REF] Cobano | Distributed thermal identification and exploitation for multiple soaring UAVs[END_REF], where data from multiple autonomous soaring aircraft is combined in order to gain information on remote updraft positions. While this technique sounds promising at first glance, an updraft might have vanished once the glider that received the information arrives. To this date, no literature can be found on active remote updraft estimation, i.e updraft estimation that is conducted on-board a UAV.

Vision-based remote updraft estimation, as studied in this thesis, deals with the estimation of the updraft position, i.e. the position of the updraft-capping cloud, and also the updraft strength and vanish time. This is a computer vision/image processing and non-linear estimation problem. Both these fields are discussed in the following.

Image Processing

While the focus of this thesis does not lie on image processing, a Master thesis project was conducted in parallel for this purpose [START_REF] Halioui | Conception d'algorithme de traitement d'image pour navigation autonome d'un drone planeur[END_REF]. For the sake of completeness, some important contributions related to the identification of clouds in images are yet summarized.

Applying computer vision techniques for cloud identification and classification is a challenging task. Clouds generally differ in size, color, texture, and geometry. Ground based or near ground based imagery, i.e. images that are taken from a perspective where the camera is located underneath the cloud contain blue pixels which originate from the sky as well as "white" pixels originating from clouds. In order to identify geometrical features of a cloud, cloud pixels have thus to be separated from non-cloud pixels. Then, cloud contours/geometrical properties can be identified.

Unlike the blue sky, clouds scatter the blue and red visible light more equally. Therefore, red/blue ratio filters have been employed for separating cloud from blue sky pixels. This approach was followed by the authors of [START_REF] Chow | Intra-hour forecasting with a total sky imager at the UC San Diego solar energy testbed[END_REF][START_REF] Johnson | A multistation set of whole sky imagers and a preliminary assessment of the emerging data base[END_REF] with a constant filter threshold. Others [START_REF] Koehler | Status of the whole sky imager database[END_REF] adopt this technique with a varying filter threshold, thus adapting the filter performance to different lightnings. More recently, cloud pixel separation techniques based on the Hue Saturation Value (HSV) color model have also been suggested. Compared to the RGB model, it provides more robustness to changes in the lightning [START_REF] Radovan | Predictions of cloud movements and the sun cover duration[END_REF]. The authors of [START_REF] Radovan | Predictions of cloud movements and the sun cover duration[END_REF] provide a method for performing cloud contour detection making use of the Suzuki Abe algorithm, as introduced in [START_REF] Suzuki | Topological structural analysis of digitized binary images by border following[END_REF].

Vision-Based Target Localization

Given vision, position, and attitude measurements, the cloud position and its geometric properties can be estimated. The position estimation process with a single camera is a monocular vision-based navigation problem. The idea is to take imagery data of the target from different vehicle locations and attitudes. Thereby, the target location can be computed by triangulation. This target localization problem, also known as bearings-only target localization, is wellknown and has its origins in the maritime sector. More specifically, in sonar applications for submarines.

Despite the prominence of this problem, the function that relates measurements to the target position is highly non-linear, which can be challenging when dealing with observer-trajectories of low observability.

The Extended Kalman Filter (EKF) is a widely used filter to cope with these nonlinearities. The author of [START_REF] Ponda | Trajectory Optimization for Target Localization Using Small Unmanned Aerial Vehicles[END_REF] for example relies on the EKF for estimating the position of stationary targets.

Others [START_REF] Watanabe | Vision-based obstacle avoidance for UAVs[END_REF] adopt the EKF for vision-based estimation of moving obstacles occuring in the flight path of a UAV. However, the EKF is sensitive to initialization and measurement errors and prone to both covariance collapse and other filter instabilities, as was reported in [START_REF] Aidala | Kalman filter behavior in bearings-only tracking applications[END_REF].

For reducing the fundamental drawbacks of the EKF, researchers [START_REF] Der Merwe | The Unscented Kalman Filter for nonlinear estimation[END_REF] have suggested the socalled Unscented Kalman Filter (UKF) where the covariance is propagated through sampling, circumventing errors introduced by linearization.

Others [START_REF] Ponda | Trajectory Optimization for Target Localization Using Small Unmanned Aerial Vehicles[END_REF] have focused on the design of trajectories that maximize the observability of the target, thus mitigating the instability of the EKF.

Yet other researchers [START_REF] Civera | Inverse depth parametrization for monocular SLAM[END_REF] have focused on increasing the measurement equation linearity by using inverse depth parametrization, thereby also minimizing the likelihood of filter collapses during initialization with large depth uncertainties.

Flight Control Law Design

Control laws for autonomous thermal soaring generally deal with the calculation of bank angle inputs that maximize the climb rate of the UAV glider, when circling within an updraft. In [START_REF] Andersson | Cooperating UAVs using thermal lift to extend endurance[END_REF], commanded bank angles are calculated by taking the vehicle's specific energy rate and its derivative as control feedback signal. In [START_REF] Andersson | Stability of a thermal centering controller[END_REF], a theoretical stability analysis of this controller is provided. Flight test results of the same controller are presented in [START_REF] Andersson | Autonomous soaring: Flight test results of a thermal centering controller[END_REF]. A slightly different control approach is adopted by the authors of [2], where the bank angle command is not only function of the vehicle's energy rate, but also of the estimated updraft center.

Concerning the glider's longitudinal motion, airspeed controllers are required to damp both the short period and the phugoid mode. The authors of [START_REF] Kahveci | Adaptive LQ control with antiwindup augmentation to optimize UAV performance in autonomous soaring applications[END_REF] apply robust adaptive longitudinal control for a soaring UAV in the presence of large model uncertainties and actuator saturations.

Guidance and Path Planning

Path planning and guidance for autonomous cross-country soaring deals with the calculation of airspeed settings and the selection of the updrafts to approach in order to reach the target.

The probably most famous publication concerning the airspeed calculation in cross-country soaring was written by physicist Paul B. MacCready, Jr., author of [START_REF] Maccready | Optimum airspeed selector[END_REF]. In this publication, an optimal airspeed setting is suggested for approaching an updraft and regaining altitude in minimum time. The airspeed setting is calculated using the glide ratio of the vehicle and the expected vertical velocity of a designated updraft. A consequence of this research is the invention of the speed-to-fly ring, a basic equipment in contemporary glider planes. Other researchers have extended this idea [START_REF] Almgren | Optimal soaring via Hamilton-Jacobi-Bellman equations[END_REF] to address the problem of optimal airspeed calculations under uncertain atmospheric conditions.

Most generally, the goal of an autonomous cross country soaring flight is to reach a remote target in minimum time with the target lying beyond the initial gliding range. In [START_REF] Kahveci | Optimal static soaring of UAVs using vehicle routing with time windows[END_REF], this problem is solved by using classical shortest path algorithms such as the Floyd-Warshall Algorithm. The authors of [START_REF] Kahveci | A heuristic search algorithm for maneuvering of UAVs across dense thermal areas[END_REF] provide a heuristic based planning solution for large updraft maps. In both these contributions, the authors assume the availability of an updraft map that contains perfect position, climb rate, and vanish time information over the entire planning horizon. The few contributions that concentrate on path planning also underline the open problem of remote updraft estimation which is vital for path planning.

Summary of Contributions

The main contributions of this work are summarized below:

Framework for Vision-Based Autonomous Cross-Country Soaring

The first contribution of this thesis is the establishment of a vision-based flight guidance and navigation system for autonomous cross-country soaring. The concept of using RGB vision for autonomous crosscountry soaring has not been studied prior to this work, and thereby provides a novelty on the field of autonomous soaring.

Remote Updraft Estimator

The second contribution of this work is a vision-based estimator for the remote estimation of sub-cumulus updrafts. The estimator not only provides the position, but also the vertical wind velocity and vanish time of updrafts. It is furthermore shown, that the estimator provides consistent state estimates which allows for an integration of estimation uncertainties in the decision-making process.

Robust Path Planner As a third contribution, this work establishes a robust path planner for vision-based autonomous cross-country soaring which takes into account the estimation uncertainty of an imperfect updraft map. The path planner is shown to be implementable on a low-cost glider UAV, using Commercials Off-The-Shelf (COTS) electronic devices.

Simulation Evaluation

For illustrating the benefits of the suggested methods, an extensive Monte Carlo (MC) simulation is conducted in a 6 DOF flight simulator with realistic glider dynamics. Performance metrics are established in order to allow a comparative study with the state of the art cross-country soaring system.

Reader's Guide

The remainder of this thesis is organized as follows:

Chapter 2 describes the design of the remote updraft estimator. This includes the description of the atmospheric and sensor models as well as the filter formulation. Monte Carlo (MC) simulations are conducted to verify the functionality of the proposed estimator in terms of convergence and consistency.

Chapter 3 explains the development of a robust path planning method for map-based autonomous cross country soaring. Therefore, a simple glider model is introduced, and a guidance law for airspeed selection is developed. The mathematical formulation of the planner is explained and simulation results are shown to illustrate the functionality of the path planner.

In Chapter 4, computer simulation results of the combined vision-based flight guidance and navigation system are presented and compared to prior systems.

Chapter 5 summarizes the results of this thesis and provides an outlook on future work.

Chapter 2

Vision-Based Updraft Map

Generation

This chapter presents a method for the vision-based, on-board, and real-time estimation of sub-cumulus updrafts. This includes the estimation of the updraft position, velocity, and vanish time, and their respective uncertainties.

Among others, the estimator relies on atmospheric and camera models, and a description of the observer dynamics. In this work, these models are selected to be as simple as possible while capturing the most significant aspects of the estimation problem. The focus in this dissertation is on general features of the remote updraft estimation and not on the particular choice of models. The estimator is designed to be flexible in the way that more complex models can easily be integrated while it remains capable of providing consistent state estimates.

This chapter is structured as follows:

In Section 2.1, the estimation problem is introduced. Section 2.2 presents an atmospheric model for updrafts and clouds which is required for the estimator design. Section 2.3 introduces the sensor models. Section 2.4 presents the estimator design for the updraft map generation. In Section 2.5, the performance of the updraft map generation algorithm is assessed and discussed by showing computer simulation results. Section 2.6 provides a brief summary of the estimator design.

Problem Statement

Consider a glider plane, deployed in an area with thermal updrafts that are capped by cumulus clouds, as illustrated in Figure 2.1. The glider is equipped with a forward pointing camera and image processing capacity, enabling a detection of cumulus clouds in the image. Here, the camera is assumed to be mounted at the vehicle's center of gravity. As a consequence, the position of the camera coincides with the position of the glider which is denoted by p. The camera's field of view is assumed not to be occluded by parts of the glider's airframe. If not explicitly indicated, all position vectors are expressed in the inertial frame. As discussed by the authors of [START_REF] Reichmann | Cross-Country Soaring[END_REF][START_REF] Pagen | Understanding the sky[END_REF], certain features of a cumulus cloud contain information on the updraft strength and the life cycle of the updraft. In this context, the cloud size plays a major role. More specifically, there is a positive correlation between the cloud size and the updraft strength [START_REF] Reichmann | Cross-Country Soaring[END_REF]. In addition, the evolution of the cloud's size indicates whether the updraft is still active or not. This is because a cumulus only grows while it is fed by an updraft and decays otherwise.

The updraft strength and vanish time, as well as the updraft position are crucial values for the decision-making in cross-country soaring as they enable the pilot/autopilot to evaluate the mission-relevant energy that can be gained from the updraft. Rather than randomly hitting updrafts, the autopilot can then generate more efficient flight trajectories and guide the vehicle to the target with less risk of mission failure [3].

In the following, the state of an updraft is denoted x u and includes the seeked parameters, i.e. its position p u , its vanish time t v , and its velocity w u,rz

x u =      p u t v w u,rz      , (2.1)
where the updraft position is noted p u = p ux p uy p uz T . The term w u,rz specifies the average updraft velocity of an updraft column with bell-shaped velocity distribution where the average refers to both the updraft radius r and the altitude z.

Given the individual updraft states x u,i , the state of an updraft map X u is defined to contain all n perceived updrafts such that

X u =         x u,1
x u,2 . . .

x u,n         . (2.
2)

The generation of such a map is an unexplored research subject and a complicated task. This is for various reasons.

For a small and low-cost realization of an estimator, vision is the only sensor that can be employed to actively identify and estimate remote updrafts. A Light Detection And Ranging (LI-DAR) based estimation would conflict with the size, weight, and cost limitations. To determine the position of cumulus clouds, UAV state estimates are combined with information gathered by vision sensors. This nonlinear, bearings-only estimation problem is sensitive to the relative motion between the observed cloud and the vehicle. Additionally, in autonomous crosscountry soaring, the flight trajectories of the glider are chosen to optimize the cross-country soaring velocity, but not the estimation observability. This fact complicates the position estimation since clouds might be poorly observable. Furthermore, the low-cost realization of such a system cannot rely on a 360 degree camera and thus the camera features a limited field of view.

Hence, vision measurements for tracked clouds are not constantly available which degrades the estimation performance. Additionally, the benefits of stereo vision can hardly be exploited due to the small ratio between the vehicle's size, i.e. the baseline, and the distance to the clouds which typically is around 3-5 km. Vision measurements and an estimate of the cloud's position can be transformed into the current diameter of the cloud. Observing the evolution of this parameter allows to gather information on updraft strength and lifespan [START_REF] Reichmann | Cross-Country Soaring[END_REF]. Imperfections of the image processing algorithms however not only degrade the estimate, but also increase the initialization time which reduces the timespan available for the UAV to benefit from an updraft.

Atmospheric Model

On days with favorable meteorological soaring conditions, updrafts create so-called cumulus clouds, when rising parcels of warm and humid air reach the Cumulus Condensation Level (CCL). As discussed in [START_REF] Reichmann | Cross-Country Soaring[END_REF], the clouds persist after the thermal has dissipated, and vanish progressively. Thus, the presence of a cumulus not necessarily indicates a thermal beneath its base. Therefore, it is important for human glider pilots to study features of clouds to differentiate clouds that can be beneficial from those that are not.

That being said, the design of a vision-based remote updraft estimator relies on an atmospheric model which describes the relation between the updraft state and cloud characteristics.

However, to this date, there does not exist a simple and suited model for this purpose in the literature. Therefore, in this work, a simple model is derived which makes the design of an estimator possible by capturing the most significant relations between the updraft state and clouds. The design of more precise atmospheric models that still balance complexity and performance is identified as future work and will be discussed in Chapter 5.

The purpose of the following section is thus the description of a model which relates the states of an updraft x u to the state of its related cloud x c

x u = f (x c ).

(2.3)

Thermal Updrafts

Convective mixing of the planetary boundary layer results in vertically oriented wind columns, referred to as thermal updrafts [START_REF] Wallace | [END_REF]. In this work, we rely on the model from [START_REF] Allen | Updraft model for development of autonomous soaring uninhabited air vehicles[END_REF] to describe the vertical velocity profile of an updraft.

Velocity profile:

According to Allen's updraft model [START_REF] Allen | Updraft model for development of autonomous soaring uninhabited air vehicles[END_REF], the velocity profile of an updraft at a given altitude features a bell shape, i.e. the updraft velocity has a peak at the core and decreases with the radius. At a given altitude z, the averaged updraft velocity of all updrafts wu,r depends on the altitude z according to wu,r (z) = z z i

1 3 1 -1.1 z z i w . (2.4)
where w stands for the vertical velocity scale and z i is the convective mixing layer thickness.

Notice that the index r specifies the average over the radius, i.e. from -r to r, and the bar stands for the average of all updrafts.

Average updraft velocity over altitude: Building upon Equation (2.4), the average updraft velocity wu,rz is introduced as the average of the velocity wu,r over a given altitude range,

i.e. between some minimum altitude z min and the cloud base altitude z cb . It is obtained by integrating Equation (2.4) over the altitude z and dividing by the altitude range z cb -z min wu,rz =

z cb z min wu,r (z)dz z cb -z min = w z cb -z min   3 4 z 4 3 z 1 3 i - 33 70
z 7 3 z 4 3 i   z cb z min . (2.5)
Note that in this notation the index rz means averaged over both the radius and the altitude. To model variations for individual updrafts, we can scale wu,rz with a perturbation gain c w u,rz,i = c i wu,rz .

(2.6)

Vanish time: As discussed by the authors of [START_REF] Kahveci | Robust Adaptive Control For Unmanned Aerial Vehicles[END_REF], the updraft velocity of a thermal barely varies with time. However, updrafts feature a limited duration. In order to model these characteristics, the perturbation gain c is defined to be constant at some positive value during the active phase of the thermal, and drops to zero at the time instance t v where the thermal van-

ishes c =      constant > 0 , if t 0 ≤ t < t v 0 , otherwise.
(2.7)

Position:

The position of an updraft is defined as the updraft's horizontal center position at the cloud base altitude z cb

p u =      p ux p uy z cb      . (2.8)
In this work, we refer to weather scenarios with no horizontal wind. Therefore, the horizontal position derivative of the updraft is null. As a consequence, updrafts can be modeled as nontilted, non-bended and torsion free.

The cloud base altitude slowly rises during the morning and descends with the nocturnal inversion [START_REF] Wallace | [END_REF]. Updrafts would thus exhibit a slow vertical movement. This phenomenon can however be neglected in the updraft model, since the duration of an updraft (∼ 20 min according to [START_REF] Allen | Updraft model for development of autonomous soaring uninhabited air vehicles[END_REF][START_REF] Pagen | Understanding the sky[END_REF]) is rather short, compared to the rise/sink time of the cloud base level.

That being said, the derivative of the updraft position vector is noted as

ṗu = 0 3×1 . (2.9)
Remark on the notation: In the remainder of this work, a matrix 0 i×j defines a null matrix of dimension i × j, and a matrix I i×j can be a pseudo identity matrix, such that

I (i, j) =      1 , if i = j 0 , otherwise. (2.10)
To give an example, I 2×3 is defined as

I 2×3 =   1 0 0 0 1 0   . (2.11)

Clouds

The state of a cloud x c can be represented by two uncorrelated vectors. These are the cloud position vector p c , and a vector d that contains parameters which describe the evolution of the cloud size w.r.t. time As shown in [START_REF] Pagen | Understanding the sky[END_REF] and further illustrated in Figure 2.2, the cloud base altitude z cb is the line intersection of the Dry Adiabatic Lapse Rate (DALR) and the dew point. According to the authors of [START_REF] Lawrence | The relationship between relative humidity and the dewpoint temperature in moist air: A simple conversion and applications[END_REF], the cloud base altitude z cb can be calculated according to

x c =   p c d   . ( 2
z cb = T 0 -T d,0 ∂T ∂z -∂T d ∂z ≈ 125 m K (T 0 -T d,0 ), (2.14)
where T and T d are the temperature of the air and the dew point respectively. The parameters T 0 and T d,0 are the respective temperatures at the ground, as shown in Figure 2.2. Only tiny variations of the temperature differences between the air and the dew point can be found when flying over regions with constant geological and orographic structure. These variations are considered to be non-significant in the mission operation zone. Thereby, the cloud base altitude z cb in the observed region can be considered as uniform for all clouds. This phenomenon is illustrated in Figure 2.2, and additionally depicted in Figure 2.3. Furthermore, this phenomenon can be exploited in the design of the cloud position estimator, as will be elaborated in Section 2.4.3. of the sky [START_REF] Pagen | Understanding the sky[END_REF]. As explained earlier, the size properties of these clouds contain information on both the updraft strength and vanish time. According to the author of [START_REF] Ackerman | Meteorology: Understanding the Atmosphere[END_REF], one characteristic only one has to be observed in order to retrieve the seeked parameters. Selecting the width rather than the height as measurement has two advantages. When observing clouds that are behind one another w.r.t. the camera, the front cloud can occlude the height of the cloud that lies behind, but not its width, when the observer is supposed to be located below the cloud base altitude. Since we want to extract the size property of the cloud from images that are recorded below the cloud base altitude, it is hence more promising to refer to the width.

Furthermore, cumulus clouds feature rather flat bases which facilitates the detection of this characteristic in an image.

As discussed by the authors of [4], the base of a cumulus cloud can be approximated by a circle with diameter d. The diameter grows while the updraft is active (t < t v ), and decays once it starts dissipating (t > t v ). Inputs to the process that describes the cloud width over time also include physical quantities such as the local humidity or air pressure. However, for a visiononly remote updraft estimation, measurements of these meteorological parameters under the remote cloud are not available. Thus, the on-board model has to be stated as a function of time only, as schematically illustrated in Figure 2.4. [4] In this work, a simple polynomial function is used for modeling the cloud size dynamics d(t) = e 1 t 2 + e 2 t + e 3 with e 1 < 0.

( (2.17)

Relation between Updraft Lifespan, Strength, and Cumulus Clouds

A decaying cloud indicates that its driving updraft has vanished. Hence, the vanish time t v corresponds to the time when the cloud's size begins to decrease, as shown in Figure 2.4. Given Equation (2.16), the vanish time is defined by the abscissa of the cloud evolution parabola's vertex

t v = - e 2 2e 1 . (2.18)
Stronger updrafts, i.e. updrafts that feature a higher updraft velocity generally result in a more significant cloud size. This phenomenon can be modeled by a linear relation between the cloud size amplitude d max and the updraft velocity

w u,i = κ d max,i , (2.19)
where the parameter d max is related to the coefficients of the parabola according to

d max = e 3 - e 2 2 4 e 1 > 0 (2.20)
and κ is some constant defined by local meteorological weather conditions and supposed to be known.

The updraft-cloud relation is then summarized as

x u = p T u t v w T = f (x c ) = p T c -e 2 2 e 1 κ e 3 -e 2 2
4 e 1 T .

(2.21)

Sensor Modeling

Camera Model

A pinhole camera model, as illustrated in Figure 2.7, is used to conduct the perspective projection from objects onto the image plane. with unit meters into a vector x ip ∈ R 2 with unit pixels in the image frame is given by

x ip = I 2×3           1 x cam z      f x 0 x 0 0 f y y 0 0 0 1      C x cam           , (2.22)
where C specifies the camera calibration matrix, x 0 and y 0 are the offsets from the upper-left hand corner to the image center in pixels and f x and f y depend on both the focal length f and conversion factors S x and S y respectively

f x = f S x , f y = f S y .
(2.23)

The conversion factors are the inverse size of a pixel in the respective axis and for a quadratic pixel, S x = S y .

Image Processing

The They contain the information on the current cloud size, given an estimate of the cloud position p c . Recent advances in image processing techniques [START_REF] Radovan | Predictions of cloud movements and the sun cover duration[END_REF] demonstrate the feasibility of detecting cumulus cloud contours in near ground imagery. Therefore, in parallel to this PhD research, a master thesis project [START_REF] Halioui | Conception d'algorithme de traitement d'image pour navigation autonome d'un drone planeur[END_REF] was launched, aiming to study image processing methods for retrieving the required measurements. In order to design the image processing algorithms, a low-cost RGB camera was mounted on a remotely piloted airplane, and video data was collected on a day with fair weather cumulus soaring conditions, creating a rich data set of images.

The crucial image processing steps that resulted from the master thesis are: For identified clouds, contour detection is applied in order to identify the polygon describing the contour of the cloud, as shown in Figure 2.9d. Finally, the cloud baseline is obtained by calculating the intersection points of the cloud contours with a line that is parallel to the horizon (Figure 2.

9e).

In what follows, an image processing algorithm that provides these two feature points is considered to be available and a vision measurement y LR ∈ R 4 is introduced as

y LR =   y L y R   = h(x c ) + ν y LR , (2.24)
where the vector ν y LR stands for zero-mean Gaussian measurement noise

ν y LR ∼ N 0, Σ 2 y LR (2.25)
and the measurement equation h in Equation (2.24) relates the cloud states to measurements.

The diagonal matrix Σ 2 y LR in Equation (2.25) specifies the measurement noise covariance of the vision measurement.

Remark on the notation: All noise terms introduced hereinafter follow the definition of Equation (2.25) and are noted ν x and ν y respectively. Subscripts x and y distinguish between pro-cess and measurement. The related sub-subscripts, (e.g. ν y LR ) denote the corresponding signal.

Center Point Sensor

As will be discussed in Section 2.4.2, we separate the updraft position estimation from the estimation of the vanish time and the updraft strength. For the cloud position estimation, we rely on vision measurements y c ∈ R 2 of the cloud's center position p c . This measurement can be reconstructed from the vision measurement y LR such that y c = g (y LR ). The geometry between the camera and the three points p L , p R , and p c is illustrated in Figure 2.10. Note that the aircraft position vector p and the two points p L and p R span a plane, and p c not necessarily lies on that plane. Also note that Ψ stands for the heading to the respective point and recall that d is the current diameter of the cloud.

In order to obtain the vision measurement y c , we start by decomposing y LR as follows

  y L 1   = C p cam L p cam Lz = 1 p cam Lz CR cam←v (p L -p) , (2.26)
  y R 1   = C p cam R p cam Rz = 1 p cam Rz CR cam←v (p R -p) , (2.27)
where the matrix R cam←v transforms a vector from the vehicle frame (noted v) to the camera frame. The vehicle frame is defined as an inertial frame with its origin at the vehicle's position.

We can rewrite Equations (2.26) and (2.27), and solve for the positions of the two outer points in the vehicle frame as function of the bearings b

v L and b v R p v L = p L -p = p cam Lz R -1 cam←v C -1   y L 1   b v L =      r cos δ   cos Ψ L sin Ψ L   ∆p z      , (2.28)
p v R = p R -p = p cam Rz R -1 cam←v C -1   y R 1   b v R =      r cos δ   cos Ψ R sin Ψ R   ∆p z      . (2.29)
It is important to note that the bearings are transformed vision measurements, i.e. available quantities. Also note that ∆p z stands for the vertical distance between the glider and the cloud.

The two bearings b v L and b v R , i.e. the bearings from the glider to the feature points in the vehicle frame can be normalized as follows

bv L =      cos Ψ L sin Ψ L ∆pz r cos δ      = b v L ||b v Lxy || , (2.30)
bv R =      cos Ψ R sin Ψ R ∆pz r cos δ      = b v R ||b v Rxy || . (2.31)
In a similar way, we can define a normalized bearing to the center point of the cloud

bv c =      cos Ψ c sin Ψ c ∆pz r      . (2.32)
Given the bearings to the outer points, we calculate the corresponding headings Ψ L , Ψ R , and the third component of the bearing to the center point, scaled by cos δ

Ψ L = arctan bLy bLx , (2.33)
Ψ R = arctan bRy bRx , (2.34
)

∆p z r cos δ = bv Lz = bv Rz . (2.35)
Following Figure 2.10, the relation between the headings Ψ L , Ψ R and the angle δ is defined by

δ = Ψ R -Ψ L 2 .
(2.36)

Since we know the angle δ, the bearing bv c becomes function of measured quantities

bv c =      cos Ψ c sin Ψ c ∆pz r      =      cos Ψ L +Ψ R 2 sin Ψ L +Ψ R 2 cos δ • ∆pz r cos δ      . (2.37)
Given Equation (2.37), we can reconstruct the vision measurement of the center point by mapping the bearing bv c onto the image plane

  y c 1   = g (y LR ) = C bcam c bcam cz = CR cam←v bv c R cam←v (3, :) bv c , (2.38)
with R cam←v (3, :) specifying the third row of the matrix.

Estimator Design

The updraft position and parameter estimation includes nonlinearities in the measurement model. A variety of filters has been proposed to cope with such nonlinearities. A common property of nearly all these methods is the idea of providing a least squares estimate of the process' state. The standard approach for nonlinear estimation is the Extended Kalman Filter (EKF). As explained in Chapter 1, this filter suffers from one major flaw: If the linearized model used in the state propagation and measurement prediction steps is only a poor approximation of the true system dynamics, large errors can be introduced in the covariance transformation [START_REF] Bailey | Consistency of the EKF-SLAM algorithm[END_REF][START_REF] Shoudong | Convergence analysis for Extended Kalman Filter based SLAM[END_REF][START_REF] Huster | Relative position sensing by fusing monocular vision and inertial rate sensors[END_REF].

Decision-making in autonomous cross-country soaring relies on accurate updraft estimates and can be improved by additionally considering estimation uncertainties. However, in case of inconsistent state and covariance estimates, the decision-making is biased which can lead to suboptimal flight performances, or even outlandings. A common way to cope with this known flaw of the EKF is to add stabilizing noise to the state covariance after each update. This procedure is however both iterative and laborious and does not guarantee better results [START_REF] Bailey | Consistency of the EKF-SLAM algorithm[END_REF].

Particle filters, where uncertainties are transformed by propagating a "large" set of randomly generated points through the system equations, provide a powerful alternative to the EKF and can also cope with non-Gaussian distributions. However, as reported in [START_REF] Thrun | Probabilistic Robotics. Intelligent Robotics and Autonomous Agents series[END_REF], the real-time application of these methods is questionable, especially when dealing with states of higher dimensions. The technique of Unscented Transform (UT), as introduced by the author of [START_REF] Julier | A skewed approach to filtering[END_REF],

is yet another alternative for transforming uncertainties through nonlinear functions. The UT has similarities to the particle-based transform. The selection of particles is however not conducted randomly, and only a few number of particles are necessitated to approximate Gaussian distributions. Thereby, the UT balances between the flaws of uncertainty propagation through linearization and the required computational effort of particle filters. The authors of [START_REF] Der Merwe | The Unscented Kalman Filter for nonlinear estimation[END_REF] applied the UT concept to the Kalman Filter which led to the so-called Unscented Kalman Filter (UKF). As demonstrated in [START_REF] Langelaan | State estimation for autonomous flight in cluttered environments[END_REF], the UKF has the capacity to provide more consistent estimates than the EKF does. Reasoning from these facts, the UKF is a well-suited candidate for the remote updraft estimation, when aiming to include information on estimated updraft uncertainties in the decision-making process.

Recall on Unscented Kalman Filtering

For the sake of completeness, the Unscented Kalman Filter (UKF) algorithm, as stated by the author of [START_REF] Der Merwe | The square-root unscented Kalman filter for state and parameter estimation[END_REF], is shortly recalled in the following.

System equations: Consider a discrete process model of type

x k = f (x k-1 , u k , ν x,k ) with ν x,k ∼ N (0, Σ 2 x,k ), (2.39)
with the state vector x, the state transition function f , the system input u, and process noise ν x . The measurement model relates the state x to the output y such that

y k = h(x k , u k , ν y,k ) with ν y,k ∼ N (0, Σ 2 y,k ). (2.40)
As any Kalman filter, the UKF follows the classical prediction-correction procedure.

Prediction:

The UKF relies on a set X of (2L + 1) sigma points which is calculated at each execution

X k-1 = xk-1 xk-1 + √ L + λ P x k-1 xk-1 - √ L + λ P x k-1 , (2.41)
where P x is the estimated state covariance, and L indicates state dimension. Also note that each column of the matrix X represents one sigma point. The constant parameter λ is function of the spread of the sigma points α and the state dimension L λ = L(α 2 -1).

(2.42)

All sigma points are transformed through the nonlinear system function f

X (i) k|k-1 = f X (i) k-1 , u k , 0 for i = 1, ..., 2L + 1, (2.43)
where the superscript i denotes the i th column of the matrix and X k|k-1 is a matrix of predicted sigma points. The subscript k|k -1 in Equation (2.43) stands for a prediction at time step k, based on inputs u, states, and measurement prior to k. The predicted average state xk|k-1 and its corresponding covariance P x k|k-1 are captured using the constant weight vectors

η c i , η m i xk|k-1 = 2L+1 i=1 η m i X (i) k|k-1
(2.44)

P x k|k-1 = Q + 2L+1 i=1 η c i X (i) k|k-1 -xk|k-1 X (i) k|k-1 -xk|k-1 T , (2.45)
where Q = Σ 2 x specifies the process noise covariance. The two weight vectors are defined as

η c i =      λ L+λ + 1 -α 2 + β , if i = 1, 1 2 (L+λ)
, otherwise.

(2.46)

η m i =      λ L+λ , if i = 1, η c i , otherwise.
(2.47) with β being the secondary scaling constant. The sigma points X (i) k|k-1 are processed through the nonlinear measurement equation which yields the predicted observations

Y (i) k|k-1 = h X (i) k|k-1 , u k , 0 , for i = 1, ..., 2L + 1.
(2.48)

The mean of the predicted measurement is defined as

ȳk|k-1 = 2L+1 i=1 η m i Y (i) k|k-1 . (2.49)
Summing the measurement noise covariance R and the covariance of the transformed state, the predicted measurement covariance is

P ȳ k|k-1 = R + 2L+1 i=1 η c i Y (i) k|k-1 -ȳk|k-1 Y (i) k|k-1 -ȳk|k-1 T . (2.50)
The prediction step is completed with the computation of the cross covariance matrix

P xy k|k-1 = 2L+1 i=1 η c i X (i) k|k-1 -xk|k-1 Y (i) k|k-1 -ȳk|k-1 T . ( 2

.51)

Correction: The calculation of the Kalman gain and the update of both estimated state and state covariance are performed as follows

K k = P x k|k-1 P ȳ k|k-1 -1 , (2.52
)

xk = xk|k-1 + K k y k -ȳk|k-1 , (2.53)
P x k = P x k|k-1 -K k P ȳ k|k-1 K T k . (2.54)
Note that here, the caret on the state helps distinguishing between propagated and updated state estimate.

Estimator Structure

The functionality of the updraft map generation is to estimate the states of the clouds with a UKF, and to transform them through the updraft-cloud relation Equation (2.21) to obtain the updraft state. However, Kalman Filters require a state initialization and are prone to initialization errors.

The cloud position p c can be initialized using only one vision measurement y LR , when assuming that some initial knowledge on the cloud base altitude is available. Initial guesses on the three parabola parameters are also required. The initialization of the parabola parameters necessitates information on the sign of the cloud growth, i.e whether the cloud is in a growing, or in a decaying stage. Contrary to the cloud position state, this information cannot be deduced from one single measurement, but requires a series of consecutive measurements. To cope with these varying initialization requirements, the cloud state estimation is divided in two filters, as illustrated in Figure 2.11. The first one provides an estimate of the cloud position, and the second one estimates the cloud size parameters, i.e. the parabola. Supposing the parameters z i , w and κ are available from a meteorological weather forecast, the vectors d can be transformed into both the vanish time and the updraft strength.

The functionality of the two estimators is described in the following. As a consequence, vision measurements for individual clouds are not constantly available since the camera features a limited field of view. That is, with irregularly available vision measurements, the estimated states of single clouds cannot be updated regularly. This inherently in-creases the filter's convergence time.

In order to compensate for these unfavorable measurement dynamics effects, induced by the glider's flight trajectory, it is important to exploit all knowledge that is available on the updrafts. Therefore, instead of running an independent filter for the estimation of each cloud position and neglecting the correlation between the individual cloud base altitudes, a combined state formulation allows a simple exploitation of the fact that all clouds lie on the same altitude level z cb . It becomes thus possible to update the estimates of all cloud positions, as soon as a single vision measurement y LR is available for one of the clouds. In other terms, the position estimate of a cloud can be updated without a vision measurement for that particular cloud. That said, the combined cloud position state estimate is defined as follows

X p =            p c xy,1 p c xy,2 . . . p cxy,n z cb            , (2.55)
where the subscript p in the concatenated state vector X p refers to the position, p c xy,i is the horizontal position of a single cloud i, and z cb stands for the common cloud base altitude.

Process Model

In the case of no horizontal wind, the cloud position vector remains constant between time steps. Therefore, the process model is defined as

Ẋp = f (ν x Xp ) = ν x Xp , (2.56)
where

ν x Xp ∈ R 2n+1 .

Measurement Model

y c = I 2×3 1 p cam c,z CR cam←v (p c -p) .
If the image processing algorithm detects m clouds, the true measurement vector is

Y c =         y c,1 y c,2
. . .

y c,m         . (2.57)
Data association Data association is required to identify whether an incoming vision measurement belongs to one of the previously initialized cloud states, or originates from a newly detected cloud. A gated nearest neighbor approach, based on the Mahalabonis distance [START_REF] Mahalanobis | On the generalized distance in statistics[END_REF] is used to tackle the data association problem for vision measurements. The underlying idea is to compare the probability that a predicted measurement corresponds to an incoming measurement. This technique has proven to work reliably [START_REF] Montemerlo | FastSLAM: A Scalable Method For The Simultaneous Localization and Mapping Problem in Robotics[END_REF][START_REF] Langelaan | State estimation for autonomous flight in cluttered environments[END_REF], provided that the uncertainty of the predicted measurement is sufficiently small.

At each filter execution, a score s is calculated for the n×m constellations between m incoming measurements and n initialized cloud position states

s i,j = y c,i -ȳc,j T P ȳ-1 c,j y c,i -ȳc,j . (2.58)
Scores s that feature a lower value than a threshold s are considered as a match. Non assigned measurements are considered to result from newly detected clouds and used to initialize a new cloud position p cxy which is appended to the combined state.

State and State Covariance Initialization

p cxy = f (y c , z cb ) = I 2×3 p + z cb p z b v cz b v c . (2.59)
For the very first cloud, this procedure requires some initial estimate of the cloud base altitude z cb and its uncertainty. The initial cloud base altitude can be obtained from the weather forecast, or, if temperature and dew point measurements are available, by means of Equation (2.14). Subsequent cloud position estimates are initialized using the estimated cloud base altitude. The related state covariance matrix is obtained by applying a UT. Therefore, the 3D measurement uncertainty of the vector y T c z cb T is transformed through Equation (2.59).

Cloud Size Dynamics Estimation

For reasons discussed in Section 2.4.2, the remote updraft estimation is subdivided in two filters, the first dedicated to the estimation of cloud positions, and the second for the estimation of the cloud size parameters. In the following, the second estimator is derived for the three parameters that define the cloud's parabolic diameter evolution function, as introduced in Equation (2.16). Given an estimate of these parameters, the vanish time and updraft strength belonging to the cloud can be obtained according to Equation (2.21).

Cloud Diameter Sensor

The crucial measurement for the estimation of the cloud size vector

y d = e 1 t 2 + e 2 t + e 3 .

State and State Covariance Initialization

Kalman Filters are sensitive to initialization parameters in the way that the initial covariance and estimation error have to be consistent. In order to initialize the three parameters of the parabola, at least three constraints about the parabola have to be known a priori. A rigorous initialization of the state requires some guess on the slope of the diameter at the time the filter is initialized. This information cannot be obtained from only one diameter measurement. Such initialization problems also occur in other domains, for example radar-based target tracking, where the velocity of the target object cannot be initialized given only one distance-measurement [START_REF] Buhren | Initialization procedure for radar target tracking without object movement constraints[END_REF].

Different methods for obtaining the slope of the diameter can be thought of. One approach is to batch-fit a series of collected cloud measurements, as suggested in [START_REF] Buhren | Initialization procedure for radar target tracking without object movement constraints[END_REF]. This solution however requires an increased memory consumption. As an alternative, a recursive initialization period is proposed, where the slope of the cloud diameter is estimated by a linear Kalman Filter. The state of the linear filter is defined as x s and contains the slope µ and the cloud diameter at the beginning of the estimation d 0

x s =   µ d 0   , (2.63)
Cloud diameter measurements are predicted according to

y d = d 0 + µ(t -t 0 ), (2.64)
where t 0 is the time at the beginning of the slope estimation.

A UKF with the state Equation (2.62) is triggered on the condition that the slope µ is positive, and its covariance inferior to a threshold value γ . Given the estimated slope µ, the current diameter d, and time t as well as some guess about the maximum diameter of the cloud d max , the three parameters of the parabola are initialized according to

e 1 = - µ 4(d max -d) , (2.65)
e 2 = µ -2e 1 t, (2.66)
e 3 = d -e 1 t 2 -e 2 t. (2.67)
The initial state covariance P x d is obtained by applying an UT of the parameters d, µ and d max and their corresponding covariances through Equations (2.65) to (2.67).

Data Flow

The flowchart illustrated in Figure 2. [START_REF] Wharington | Autonomous control of soaring aircraft by reinforcement learning[END_REF]

Y LR,k Y c,k Y unin. d,k Y in. d,k X new d,k Y unin. c,k Θ k X new

Simulation Results

This section presents computer simulation results. The main purpose of the simulations is to study the filter performance in terms of estimation consistency and convergence.

Context

Kalman Filters not only provide an estimate of the state, but also the estimated state covariance.

The state covariance constitutes a metric of the accuracy of the estimation. As will be discussed in Chapter 3, it is beneficial to include this supplemental information in the path planning task, to render the decision-making process more robust. However, this requires consistent state covariance estimates. An inconsistent ("underconfident" or "overconfident") estimate can result in less efficient or high-risk trajectories. Applied to path planning for autonomous crosscountry soaring, an inconsistent estimation would result in lower cross-country speeds, or, even worse, in an increased number of outlandings. Therefore, it is important to not only study the convergence of the filter, but also its consistency.

Parameter Values

The parameters required for the estimator are the ones of the camera, and the filter/sensor settings.

The camera characteristics that were defined in Section 2.3 are specified in Table 2 The filter parameters include the process and measurement noise covariances, as well as the sampling time. Both the process and measurement noise covariance are defined, such that

Q = E ν x ν T x = Σ 2 x , R = E ν y ν T y = Σ 2 y .
(2.68)

Recall from Equation (2.25), that the diagonal matrices Σ specify the standard deviations. The corresponding parameter values are listed in Tables 2.2 and 2.3.

Parameter

Value Unit

Σ xp cxy 0 2×2 [m/s] Σ xz cb 0 [m/s] Σ xs I 2×2 0.005 0.1 T [m/s 2 ], [m/s] Σ x d I 3×3 0 0 0.02 T [m/s 3 ], [m/s 2 ], [m/s]

Evaluation Metrics

A criterion that is dedicated to study the consistency of an estimator is the so-called Normalised Estimation Error Squared (NEES), noted . Following the definition from the author of [START_REF] Yaakov Bar-Shalom | Estimation with Applications To Tracking and Navigation[END_REF], the NEES at sample k is defined as

k = xT k P -1 k xk , (2.69)
where x specifies the error between the true and the estimated state, i.e. xx. Under linear and Gaussian assumptions, the average of the squared norm of the estimation error satisfies the following hypothesis

E xT k P -1 k xk = L, (2.70)
where L is defined as the state dimension. Equation (2.70) is true in case the estimator correctly estimates the state covariance P.

The average NEES can be calculated from MC simulations. For a total of N simulation runs, the N -run average NEES is defined as

¯ k = 1 N N i=1 i k , (2.71)
where i indicates the run.

The calculated NEES can be verified by applying a χ 2 -test. The test is passed, if the NEES remains within an acceptance interval

¯ k ∈ [r 1 , r 2 ] , (2.72)
where r 1 and r 2 specify the bounds of the acceptance interval. These bounds depend on the selected significance level α such that

1 - α 2 = F χ 2 (N ×L) (r 1 ) ⇒ r 1 = F -1 χ 2 (N ×L) 1 - α 2 (2.73) α 2 = F χ 2 (N ×L) (r 2 ) ⇒ r 2 = F -1 χ 2 (N ×L) α 2 , (2.74)
where F stands for the Cumulative Distribution Function (CDF) which is defined as

F χ 2 (N ×L) (x) = x 0 f χ 2 (N ×L) (y) dy. (2.75)
While the NEES is suitable to study the consistency, the Root Mean Square Error (RMSE) provides a metric to evaluate the estimation performance in terms of convergence. The RMSE at sample k over N simulation runs is defined such that

ξ k = 1 N N i=1 xT k xk .
(2.76)

Monte Carlo Simulation Scenario

For obtaining representative estimation conditions, the MC scenario is specified by a typical thermal centering flight trajectory, as illustrated in Due to the spiral-shaped flight trajectory, clouds repeatedly appear and disappear on the image sensor, i.e. vision measurements are not constantly available. This creates realistic observation conditions in terms of measurement regularity.

It is further worth mentioning, that the observation trajectory stays representative, beyond the fact that horizontal wind is not considered. This is because horizontal wind does not change the relative motion between the glider and the observed clouds, since both updraft and clouds drift downwind with the same speed. Thereby, the relative motion is not impacted.

For the Monte Carlo simulation, initial states are imposed a priori. This allows to generate samples of the true states, according to both the initial state and state covariance, i.e. x 0 and P 0 respectively. At each run, the true state is sampled from this normal distribution which generates N different scenarios with

x ∼ N x 0 , Σ 2

x 0 P 0 = E x 0 x T 0 = Σ 2 x 0 .
(2.77) This threshold is usually reached when the cloud position estimation has converged. In other terms, the vanish time and updraft strength estimation requires the cloud position state to be converged. To account for this, two dedicated Monte Carlo simulations are conducted.

In the first one, the cloud position estimation is evaluated. For the second MC simulation run, the initial uncertainties in the cloud position states are selected from a normal distribution with less uncertainty, as noted in Table 2.6. This simulates that the position estimate has converged and the vanish time and strength estimation has been triggered. The remaining settings are the same as in the first run.

Parameter Value Unit

Σ pc xy,i 5 I 2×2 [m] Σ z cb 5 [m]
Table 2.6: Initialization parameters for the second MC run To illustrate this effect, the simulation scenario from Figure 2.14 is considered and repeated twice. In the first run, only the cloud number one is estimated, i.e. estimation updates depend on measurements of only one cloud. In the second run, all four clouds are tracked. Further simulations demonstrate the benefit of the combined cloud position state.

Consistency and Convergence

corporating uncertainties, it also provides a local planning solution in case the destination lies outside the range of the estimated map. Finally, it is shown that the suggested path planner is implementable on a low-cost, miniaturized computer platform, thus being a solution that can be embedded on-board a small glider UAV.

This chapter is organized as follows:

In Section 3.1, the map-based cross-country soaring problem is stated. Section 3.2 serves to explain the glider motion model, necessary for the design of the path planning algorithm and suggests a guidance law for the final glide. In Section 3.3, the impact of uncertain updraft estimates in cross-country soaring is explained. Section 3.4 provides an overview of planning in uncertain environments and recalls the fundamentals of Markov Decision Processes (MDPs) and dynamic programming. In Section 3.5, the cross-country soaring problem is formalized as a MDP and it is shown how its solution can be calculated. Section 3.6 presents computer simulation results that compare the proposed method to former approaches and a study of the runtime. The content of this chapter is summarized in Section 3.7.

The Path Planning Problem

As illustrated in Figure 3 This minimum altitude could for instance be defined by airspace restrictions.

The success and performance of a cross-country soaring flight mostly depend on three factors.

1. The vehicle's aerodynamic performance. High performance sailplanes with a high glide ratio can transform altitude into traveled distance more efficiently, i.e. faster and at a lower sink rate.

2. The meteorological conditions. Both the density of updrafts and updraft velocities directly impact the flight time.

3. The decision-making of the pilot/autopilot w.r.t both the flight path and airspeed. If the vehicle is not equipped with any internal propelling system, inefficient or wrong decisions not only increase the flight time, but can lead to so-called outlandings where the vehicle has to abort the flight by landing on a site which is not the targeted one.

In this thesis, we seek to improve the cross-country soaring performance by advanced decisionmaking. Therefore, we have to process information on the vehicle performance and the updraft map in order to calculate a sequence of consecutive control inputs/actions that bring the vehicle to the target in minimum time, if possible, or to decide to abort the mission and land the vehicle safely. To fulfill this task, models that describe the vehicle dynamics, and its interaction with meteorological phenomena are required.

Airspeed Selection and Glider Motion Model

For the path planning task, a motion model is required in order to propagate the glider's state.

We rely on a simplified point-mass model, for the sake of computational efficiency. As illustrated in Figure 3.1, the glider can operate in three different modes which are interthermal glide, thermal climb and final glide. In cross-country soaring, the high-level control inputs (actions) are the sequence of updrafts to visit, and the corresponding airspeed commands. For keeping the path planning problem tractable, we separate the airspeed selection from the path planning task. That is, the airspeed is not considered as a control variable in the path planning, but calculated as function of a separate guidance law. That said, this section also serves to identify suitable airspeed commands.

In the following, the vehicle's state s is defined as

s =   p t   , (3.1)
where p = [p x p y p z ] T is the glider's position and t stands for the time.

Interthermal Glide

The maneuver for harvesting an updraft's energy can be separated in two consecutive sequences. As depicted in Figure 3.1, these sequences are interthermal glide and thermal climb.

The interthermal glide is the straight and steady flight phase between some state s and the arrival state s a underneath the cloud of the targeted updraft. This sequence is illustrated in

V g -v si   , (3.2)
where V g and v si are the horizontal speed over ground and the sink rate respectively, the horizontal distance to the updraft is defined by q u = ||p uxy -p xy ||, and the term p uxy -p xy /q u is the direction vector of the glider in the horizontal plane. Assuming no wind, the relation between the airspeed V a and the ground speed V g is defined as illustrated in Figure 3.3 where γ denotes the flight path angle. The sink rate v si is defined positive downwards. For soaring V a airspeed V g ground speed γ flight path angle Figure 3.3: Glider kinematics in the vertical plane during glide planes, the sink velocity v si is small, compared to the airspeed V a . Therefore, both airspeed and ground speed approximately coincide

V g = V a cos γ ≈ V a , γ < < 1. (3.3)
In this work, the MacCready airspeed V M C , introduced by the author of [START_REF] Maccready | Optimum airspeed selector[END_REF], is selected as longitudinal guidance command during interthermal glide. It is the optimal airspeed for minimizing the cumulated time of an interthermal glide and its succeeding climb back to the initial altitude at the beginning of the glide-climb-sequence. The derivation of V M C is shortly recalled. Consider the glider-updraft configuration as introduced in Figure 3.2. Assuming the climb rate of the vehicle in the approached updraft w is known, V M C can be derived as follows:

The cumulated time t cum of a glide-climb sequence is the sum of the glide time t g and the climb time t c to the initial altitude p z

t cum = t g + t c = t g + ∆ z w = q u V a tg + q u V a v si w tc = q u V a 1 + v si w . (3.4)
with w specifying the glider's climb rate during thermal climb. The vehicle's polar curve that relates airspeed to sink rate can be stated as

v si = aV 2 a + bV a + c, (3.5)
where a, b, c are the known coefficients. We can substitute the sink rate in Equation (3.4), yielding

t g + t c = q u V a 1 + aV 2 a + bV a + c w = bq u w + q u V a w aV 2 a + c + w .
The optimal airspeed setting is therefore the solution of

V M C = min Va q u V a w aV 2 a + c + w . (3.6)
Differentiating with respect to the airspeed, and solving yields the MacCready airspeed

V M C = c + w a . (3.7)
Note that the MacCready airspeed solely depends on the climb rate and the vehicle's aerodynamic performance.

V M C = 1.21+w 0.0132 w [m/s] V M C [m/s]
s a = f (s, x u ) = s +      pu xy -pxy qu V M C -v si (V M C ) 1      q u V M C =      p uxy p z,a t a      , (3.8)
where the arrival altitude p z,a can be expressed as function of the arrival time and the Mac-Cready speed setting

p z,a = p z -v si (V M C ) (t a -t) ,
with

t a = q u V M C + t.

Thermal Climb

As further illustrated in Figure 3.2, the thermal climb phase is the transition from the arrival state s a to the exit state s e . From the updraft map, we know the average updraft velocity of an updraft w u,rz , as introduced in Equation (2.6). The climb rate w does however not only depend on the updraft velocity, but also on the glider's sink rate and its position relative to the updraft center. According to the updraft model from [START_REF] Allen | Updraft model for development of autonomous soaring uninhabited air vehicles[END_REF] the updraft velocity profile w u is represented by a bell curve. Therefore, the updraft velocity decreases from the chord line of the updraft outwardly. Assuming the glider to be perfectly centered on a circle aroud the updraft's center where the radius of the circle is small compared to the dimension of the updraft, the vehicle will primarily be impacted by the updraft's core velocity w p . From [START_REF] Allen | Updraft model for development of autonomous soaring uninhabited air vehicles[END_REF], the relation between w p and the average updraft velocity w u,r at a given altitude is

w p = 3w u,r r 3 2 -r 2 2 r 1 r 3 2 -r 3 1 , (3.9)
where r 1 and r 2 are the inner and the outer radius of the updraft. As further explained in [START_REF] Allen | Updraft model for development of autonomous soaring uninhabited air vehicles[END_REF], both radii depend on the altitude and the convective mixing layer thickness z i . However, for the sake of simplicity, the radii are assumed to remain constant over altitude. For typical convective mixing layer thicknesses of z i = 2000 m, the fraction in Equation (3.9) is of the order of 0.7. We can therefore state a simplified relation between the average updraft velocity w u,rz and the average core velocity w p,z w p,z = 3 • 0.7 w u,rz .

(3.10)

Compensating w p,z for the glider's sink rate during climb, the climb rate w of the glider is stated as

w = w p,z -v si (3.11)
Then, the position derivative in thermal climb can be defined as

ṗ =      0 0 w      . (3.12)
Note that Equation (3.12) is a simplification, serving as computationally efficient motion prediction model. In reality, the glider performs a spirally shaped climb trajectory. The altitude at the exit state can be noted as

p z,e = p z,a + ∆z,
where ∆z is the gained altitude. In this work, we consider two conditions for terminating a climb. First, when the glider reaches the cloud base altitude z cb . Second, when the updraft vanishes. That being said, the gained altitude ∆z is defined as

∆z = min(t v -t a , z cb -p z,a w )w, (3.13)
where t a stands for the arrival time at the updraft. With Equation (3.13), the exit state s e is a function of both the arrival state s a and the updraft vector

x u s e = f (s a , x u ) = s a +         0 0 w 1         ∆z w . (3.14)

Final Glide

The final glide is the transition between some state s and a goal state s G such that the horizontal position of the glider corresponds to the horizontal position of the destination point, i.e. p xy = p Dxy . Recalling Equation (3.3), the duration t f of the final glide can be stated as the ratio of the horizontal distance to the destination and the airspeed setting

t f = q D V a . (3.15)
As an optimal airspeed for the final glide, we seek to find the airspeed setting V a that minimizes the flight time t f of the final glide such that the minimum altitude constraint is not violated. Mathematically, this can be stated as

min V a,min ≤Va≤Va,max -v si (Va)t f +pz>z min t f , (3.16)
where V a,min and V a,max define the airspeed envelope. Note that for a glider, t f is minimized when arriving the lowest at the destination. Therefore, the following holds

t f = p z -z min v si = p z -z min aV 2 a + bV a + c , (3.17)
where p z is the glider's altitude at the beginning of the final glide. Rearranging, the seeked optimal airspeed setting for final glide, V a is obtained by calculating the zeros of the following quadratic function

0 = aV 2 a + b + z min -p z q D V a + c with V a = - 1 2a b + z min -p z q D + 1 2a b + z min -p z q D 2 -4ac. (3.18)
Notice that the final glide law can easily be extended to incorporate horizontal wind, as shown in [START_REF] Stolle | A vision-based flight guidance and navigation system for autonomous cross-country soaring UAVs[END_REF]. The airspeed setting V a , as function of the excess altitude p z -z min and the remaining horizontal distance to the destination q D , is illustrated in Figure 3.5 for the Cularis glider model from Appendix A.

Also note that here, V a is bounded by the airspeed of the best glide V a,min = 10 m/s and the maximum airspeed of V a,max = 25 m/s respectively. In other terms, if the destination cannot be reached without undercutting the minimum altitude, the selected airspeed for final glide is the airspeed of the best glide. Conversely, a maximum airspeed of 25 m/s is tolerated.

Given the airspeed for the final glide, we can introduce goal states s G ∈ s|p xy = p Dxy . Goal 0.5 1

•10 4 500 1,000 10 20 q D [m] p z -z min [m] V a [m/s]
s G = f (s, p D ) = s +      p Dxy -pxy q D V a -v si (V a ) 1      q D V a . (3.19)

Cross-Country Soaring in an Uncertain Environment

Description

Different flight paths have to be compared for identifying the shortest path in time to the target. Therefore, consider Figure 3.6. Suppose the pilot/autopilot has to make a decision: One option is to approach the updraft, regain altitude and then to glide to the destination. The second option is to directly glide to the destination. s current state s a arrival state s G,2 state at the destination for route 2 s e exit state Figure 3.6: The effect of p z,e on outlandings Furthermore, suppose the first option 1 appears to be faster in time, since the vehicle can approach the updraft with an increased airspeed, regains altitude and initiates a final glide with increased airspeed towards the destination. In contrast, the second option 2 is slower in time since no more potential energy is gained on the way to the destination and the vehicle has thus to fly at a lower airspeed for not arriving at the destination below the minimum altitude.

For the decision-making, the vehicle's current state s, the glider's polar curve parameters, and the estimated updraft parameters xu = pT u tv ŵ T are given, as introduced in Section 2.1.

Additionally, the uncertainty which is related to the updraft state is known in terms of the corresponding covariance matrix P xu . Given Equations (3.8), (3.14) and (3.19), the pilot/autopilot can determine the state s G and evaluate both the cumulated flight time and the arrival altitude for a given flight path. However, all updraft parameters are uncertain, rendering the predicted state s T inherently uncertain. As was shown in Chapter 2, the updraft position tends to converge significantly faster than both the estimates of the vanish time and the updraft velocity.

We therefore only consider uncertainties related to these parameters, i.e. the vanish time and updraft velocity. The uncertainties of these parameters enter the state propagation process through Equations (3.13) and (3.14). If we start propagating the vehicle's position from s, the first uncertain state will be s e . Note that the arrival state s a is deterministic since both altitude and time depend on the MacCready airspeed which is known a priori for a given updraft.

Also note, that the uncertainties are considered to remain constant over time. As horizontal position uncertainties are not considered, one can separate states s according to

s =   s d s p   , (3.20)
with

s d =   p x p y   s p =   p z t   , (3.21)
where the subscripts d and p stand for deterministic and probabilistic respectively. Since the estimation error is considered Gaussian, s p is a random vector

s p ∼ N (s p , P sp ) . (3.22)
Reconsidering the example from Figure 3.6 both the altitude and time of arrival at the destination are uncertain when deciding for route 1 . In other terms, this particular route might be faster but also bears the risk of an outlanding. The state uncertainty primarily depends on the spatial constellation between the vehicle, the updraft and the destination, as well as the estimation parameters. While the constellation can obviously not be changed, it is possible to integrate the estimation uncertainties in order to thoroughly evaluate whether a decision is beneficial or not. This explains the importance of considering the vanish time and the updraft strength uncertainties in the path planning problem for VBACCS. Integrating these uncertainties, the planning task from Section 3.1 becomes a constrained and time-dependent Stochastic Shortest Path (SSP) problem.

A Numerical Example

Suppose the following inputs are given to the simple decision-making process from

P xu =            0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 (200s) 2 0 0 0 0 0 (0.4m/s) 2            , p Dxy = 7700 m 0 m T , z min = 500 m,
where the updraft state is supposed to be Gaussian, such that x u ∼ N (x u , P xu ).

f s Gp = 1 2π|P s Gp | 1 2 e -1 2 s Gp -s Gp T P s Gp -1 s Gp -s Gp , (3.24)
with | •| standing for the determinant of a matrix. The probability of an outlanding is the result of the following integral

P (p z,D < z min ) = ∞ -∞ z min -∞ f s Gp dp z,D dt = 0.47. (3.25)
Thus, on the regarded example, the vehicle has a 47% risk of of an outlanding when deciding to approach the destination via the designated cloud. In contrast, by means of Equation (3.19), the flight time of a direct final glide is 611 s, i.e. approximately 19 s slower, however without the risk of an outlanding. MDPs are sequential decision models and were introduced by the authors of [START_REF] Bellman | A Markovian Decision Process[END_REF] and [START_REF] Howard | Dynamic Programming and Markov Processes[END_REF]. According to the authors of [START_REF] White | A survey of applications of Markov Decision Processes[END_REF], the broad application of MDPs ranges from finance, over agriculture, to inventory and production. With a trend toward more efficient computer processors, MDPs have also gained the attention of researchers in robotics. For an example, the authors of [START_REF] Ng | Autonomous inverted helicopter flight via reinforcement learning[END_REF] perform autonomous inverted helicopter flight, relying on reinforcement learning and MDPs. Others have relied on MDPs for UAV flocking [START_REF] Quintero | Flocking with fixed-wing UAVs for distributed sensing: A stochastic optimal control approach[END_REF] or target following [START_REF] Baek | Optimal path planning of a target-following fixed-wing UAV using sequential decision processes[END_REF]. A work with similarities to the VBACCS problem is provided by the authors of [START_REF] Al-Sabban | Wind-energy based path planning for unmanned aerial vehicles using Markov Decision Processes[END_REF], where the framework of MDPs is utilized for UAV energy harvesting in horizontal winds.

Recall on Markov Decision Processes and Dynamic

For systems with a large number of states and infinite planning horizons, MDPs tend to be computationally expensive. If the system's uncertainty is known in advance, the decision-making problem can be solved off-line where the calculation time is not predominant. However, for VBACCS, there is no way of quantifying the uncertainty beforehand. This is because updrafts occur randomly over the region the vehicle operates in, which is an unpredictable process.

Therefore, the decision-making problem has to be solved in real-time, and preferably onboard the vehicle. As compared to many other applications where the problem size has to be reduced, in VBACCS, both the comparably small problem size and a naturally limited horizon allow a near-optimal onboard solution calculation in short time. That said, the MDP framework is adopted for tackling the VBACCS planning problem.

In the following, a special class of MDPs, namely Stochastic Shortest Path MDPs [START_REF] Kolobov | Planning with Markov Decision Processes: An AI Perspective[END_REF] will be recalled which is particularly suited to model the VBACCS planning problem.

Stochastic Shortest Path Markov Decision Processes

This work relies on the SSPMDP definition, proposed by the author of [START_REF] Kolobov | Planning with Markov Decision Processes: An AI Perspective[END_REF]. The SSPMDP is represented by a tuple M = S, A, P, C, G,s 0 , K . The components of the tuple are defined as:

• S indicates the set of all possible discrete system states.

• A is a finite set of discrete actions that the agent can execute. A(s) ⊆ A defines the subset of applicable actions at s. That is, A(s) = {∀a ∈ A|I(s, a) is true } where I(s, a) is true when a is applicable at s and false otherwise.

• P : S × A × S stands for the stochastic transition model where P(s |s, a) is the probability of transitioning from state s to state s when executing the action a in state s.

• C : S × A × S → [0, ∞) is the cost model such that C(s |s, a) defines the cost of transitioning from s to s when executing action a in state s. C(s |s, a) is strictly positive whenever s is not a goal state.

• G ⊆ S defines the set of all goal states s G , also known as absorbing states, such that for all s G ∈ G, for all a ∈ A, and for all s / ∈ G the following holds: P (s / ∈ G|s G , a) = 0, and

C (s G |s G , a) = 0.
Reaching any of the states in G terminates the execution.

• s 0 is the initial state.

• K defines a finite set of decision epochs such that K . . = {0, . . . , N }. Elements of K specify the current decision epoch with running index k. Throughout this thesis, a convention is adopted where the last decision is made at the epoch N -1.

In comparison to the more general case of MDPs, the SSPMDP features some major particularities. It deals with costs and not with rewards. Therefore, the problem is a minimization and not a maximization problem. Additionally, SSPMDPs contain absorbing goal states, i.e. states

where actions have no effect on the cost/reward.

The MDP related Problem for Finite Horizons

The problem which is related to a SSPMDP consists of finding a policy, i.e. a sequence of actions π = (π 0 , π 1 , ...π N -1 ), that enables to bring the agent from some current state to a goal state, while minimizing a given criterion. Note that the element π k (s k ) of a policy specifies the action to be taken in node s at decision epoch k.

Let the value function V represent a metric for the utility of a policy. More specifically, V π specifies the expected total transition cost under some policy π, when starting in a state s 0 in the first decision epoch.

V π (s) = E N -1 k=0 C (s k+1 |s k , π k (s k )) s 0 = s . (3.26)
The problem of finding the optimal policy π can then be stated such that

π = min π V π (s). (3.27)
Different policies have to be evaluated and compared against each other, in order to identify the optimal policy π . However, it is not trivial to calculate the cumulated expected cost from Equation (3.26) in a forward way. Instead, recursion can be applied to determine the expected cost over N -k remaining decision epochs

V π k (s) = s ∈S P s |s, π(s) C s |s, π(s) + V π (s ) .
(3.28)

Backward Induction for Solving SSPMDPs

In the following, all values V N corresponding to all states of the set S N ∈ S, i.e. states at decision epoch N are supposed to be given. Since we seek to identify the next action to execute for a state with N -1 remaining actions to take, we can apply the recursion from Equation (3.28)

to project the values of the states S N back in time, i.e. in reversed chronological order while minimizing the expected future cost and extracting the corresponding policy. The underlying idea is the Bellman principle. It states that an optimal policy has the property that whatever the initial state and decision are, the remaining decisions must constitute an optimal policy with regard to the state resulting from the first decision. The identification of the optimal values V k leads to the corresponding optimal policy π which is defined as Note that backwards induction can also be thought of as a version of the more widely known value iteration algorithm [START_REF] Thrun | Probabilistic Robotics. Intelligent Robotics and Autonomous Agents series[END_REF], applied to finite horizon processes.

Starting in decision epoch

π k (s) = arg min

Definition of Actions

In this work, two different types of actions are considered. These are actions of type a α : transition to thermal and actions of type a β : transition to target. The first action is defined by a sequence of an interthermal glide (Section 3.2.1), followed by a thermal climb (see Section 3.2.2 for an explanation). The applicability function for actions a α is defined such that

I (s, a) |a ∈ A α =      true, if p z,a ≥ z min and |∆Ψ| ≤ 0.5π false, otherwise, (3.32)
where A α is the set of all actions a α . Recall that p z,a is the glider's altitude when arriving at the designated updraft. That is, following Equation (3.32), an updraft can be approached on condition that the glider reaches the updraft position without violating the minimum altitude constraint. Furthermore, course deviations ∆Ψ w.r.t. the target are limited to ±90 • .

Subsequently, the glider can only approach updrafts that lie in the sector between its current position and the target.

In contrast to transition-to-thermal actions, transition-to-target actions are defined as always applicable

I (s, a) |a ∈ A β = true, (3.33)
where the set A β contains all actions of type a β . Though the applicability for transit-to-target actions is true from any state, this does not mean that the glider will reach the target without undercutting the minimum altitude constraint from any state.

State Transition Model

In this work, it is explicitly distinguished between the deterministic and the probabilistic part of the state transition model. The deterministic part relates mean states to mean outcomes over actions. The probabilistic part of the transition model deals with the calculation of discrete state transition probabilities.

We recall, that the framework of MDPs can also be used to calculate a deterministic planning solution. To do so, it suffices to ignore the stochastic part of the transition model. This is done in the evaluation of this chapter, see Section 3.6, in order to illustrate the benefits of probabilistic path planning in VBACCS, compared to deterministic path planning.

Deterministic State Transition Model

The transition model for a α actions takes a state s as input and predicts an exit state s e based upon the glider motion model from Section 3.2.

s e =      p u,xy p z -v si (V M C ) qu V M C + ∆ z t + qu V M C + ∆z w      . (3.34)
In the sequel, we will define the set of goal states S G = {S|p x = p Dx , p y = p Dy }. Transition to target actions a β are then defined to relate states s to goal states s G ∈ S G according to where the integration region R is defined to not only depend on the area B of the related bin.

s G =      p Dxy p z -v si q D V a t + q D V a      . ( 3
In particular, we define the region as Thereby, if a bin is adjacent to parts of the constrained state region, all probabilities of terminating the action from the adjacent constrained region are mapped onto the bin.

R =                                            {(p z,
The composition of the bin-related probabilities is explained in Figure 3.10. The first colormap, i.e. Figure 3.10a, illustrates the probabilities, obtained when only integrating over the area. Note that the frayed border marks the outer bound of the constrained region. Figure 3.10b illustrates the bin-related probabilities following Equations (3.38) and (3.39). Consider the upper right bin for an example. According to Figure 3.10a, the probability of terminating an action from this particular bin is approximately zero, when not considering its adjacent regions.

However, following Equation (3.39), the probabilities of its adjacent regions are mapped onto it, which increases its related probability. The bin's colour shifts from white to light grey. In the same way, the probability of terminating an action from the upper center bin increases, its colour shifts from dark grey to black.

Note that there is no closed-form solution for the integral of the bivariate Gaussian distribution Equation (3.38). Rather than solving the integral iteratively, the approximation for the bivariate case, as suggested by the authors of [START_REF] Drezner | On the computation of the bivariate normal integral[END_REF], is adopted. It constitutes a computationally

efficient and yet precise solution. The discretized state with both bin-related probabilities and conditional expectations is illustrated in Figure 3.11. As expected, the bin-conditional expectations are not centered in the bin.

t
In addition, Figure 3.12 symbolically illustrates the discrete outcome states of an action and its related probabilities (from the side view).

Cost Model

A cost model C is required for defining state transition costs. In the autonomous cross-country soaring problem, the transition cost for actions a α is the sum of the glide time and the climb time to the respective altitude. Recall that this quantity can be calculated following Equation (3.4). Thus, we can state

C(s |s, a α ) = t g + t c = q u V g + ∆z w . (3.41)
The cost calculation for transitioning to the target, i.e. actions a β is illustrated in Figure 3.13 where V a stands for the constant airspeed during final glide.

s current state s G goal state γ flight path angle V a airspeed V g ground speed v si sink rate z min minimum alt.
q 0 h. dist. to outland. q r remaining h. dist. q D h. dist. to destination In case the destination can be reached without violating the minimum altitude constraint (see 1 in Figure 3.13), the cost is defined as the glide time.

In case the target cannot be reached without undercutting the minimum altitude z min (see 2 in Figure 3.13), the cost is defined as the glide time to the outlanding site plus a penalty Ψ.

For the two cases, the transition cost is summarized as

C(s |s, a β ) =      q D V a , if p z,D ≥ z min q 0 V a + Ψ , otherwise, (3.42)
where the corresponding distances are

q D = ||p Dxy -p xy ||, (3.43
)

q 0 = pz-z min -v si (V a ) V a , (3.44
)

q r = q D -q 0 . (3.45)
The penalty Ψ for an outlanding contains both a constant term k 0 and a term that depends on the remaining distance

Ψ = q r k 1 V a + k 0 . (3.46)
The dimensionless parameter k 1 serves to scale the constant final glide velocity V a . Recall that V a corresponds to the best airspeed of maximum lift-to-drag, for the case the glider cannot reach the destination. The term k 1 V a in Equation (3.46) can be thought of as the velocity of a fictitious conveyor belt that transports the vehicle from the outlanding site to the destination.

Then, the first term of Equation (3.46) stands for the transportation time it takes to bring the vehicle from the outlanding site to the destination. The second term, i.e. the time k 0 encapsulates the glider competition rule, that landing out one meter short of the target is always worse than accomplishing the mission without an outlanding [START_REF] Almgren | Optimal soaring via Hamilton-Jacobi-Bellman equations[END_REF].

Integrating the notion of an outlanding into the cost is fundamental for the purpose of this thesis. Through the penalty for an outlanding, we can balance between the risk of a mission failure and the performance in terms of flight time.

Tree Construction for Policy Optimization

Given the problem in MDP form, a decision tree featuring N + 1 layers is generated, as illustrated in Figure 3.14 where the rectangles/squares define states and the circles specify actions.

The tree is constructed in chronological order, starting from the vehicle's current state s 0 . At each node in the layers 0 ≤ k < N , the applicability function outputs transit-to-thermal actions a α for the thermals that lie in the glide range and the transit-to-target action a β . As an additional condition, the vehicle can harvest energy only once from each updraft. That is, after having

s 0,1 a β s G 1,N 1 a αr a α 1 s 2,(m×n) a β s 2,1 a β a β s G N,2 a β s G N,1 a β s G N,
k = 0 k = 1 . . . k = N -1 k = N s i,
k > 0.
In other terms, a goal state can be reached from any of the non-goal states. Recall that this does not necessarily mean that the destination is indeed reachable. However, this formalization allows to generate local planning policies that are optimal with respect to the amount of information that is currently at hand. At this point, the reader is also reminded about a key assumption in the strong definition of SSPMDPs, as provided by [START_REF] Kolobov | Planning with Markov Decision Processes: An AI Perspective[END_REF]. There exists at least one policy that will eventually terminate in a goal state, irrespectively of the initial state. With the goal state being applicable from any state, this key property of SSPMDPs is maintained.

Since the tree is constructed in a way that only transit-to-target actions can be applied from the second last layer (k = N -1), all leaf states in layer N are goal states.

The pseudo code for the tree construction is stated in Algorithm 2.

The asymptotic time complexity of the tree construction can be stated as

O |A||m × n| N -1
where the primitive operation is defined as the function which is called in line four of the algorithm.

. Planning Horizon Selection

For Vision-Based Autonomous Cross-Country Soaring (VBACCS), the planning horizon can be deduced from typical meteorological parameters and the vehicle's aerodynamic performance.

Recall from Section 3.5.5, that the vehicle can execute at most N actions to reach the goal and the last action is always a transit-to-target action. By consequence, up to N -1 updrafts can be approached.

Using vision for the generation of the map limits the perception horizon to the sensing range.

Furthermore, updrafts feature a typical lifespan of 20 min [START_REF] Allen | Updraft model for development of autonomous soaring uninhabited air vehicles[END_REF]. Thereby, the generated map is at most valid for a duration of 20 min. If we know the average time it takes to execute a transit-to-thermal action taα , we can solve for the planning horizon according to

N = 20 min taα + 1. (3.47)
As was shown in Equation (3.41), the duration of a transition to thermal action depends on the distance to the updraft, the climb rate and the vehicle's aerodynamic performance. Literature provides typical values for climb rates and the distance between updrafts [START_REF] Müller | Meteorologische Aspekte des Streckensegelfluges: Thermik von A bis Z[END_REF]. Therefore, one can determine a representative value for taα and eventually solve for N .

Simulation Results

The primary purpose of this section is to illustrate the benefits of the proposed path planner by means of a simple scenario, assuming that an imperfect updraft map is provided. An additional analysis is conducted to study the required runtime of the suggested path planner on a typical state-of-the-art computer processor.

Numerical Example

In what follows, the 20 km cross-country soaring flight illustrated in Figure 3. -2,000 0 2,000

#1 #2 #3 #4 y [m] x [m]

Updrafts

Start Destination A map with imperfect updraft parameters and related uncertainties is generated as input to the planning problem. Note that the uncertainties only refer to the vanish time and updraft climb rate, while the position is perfectly known.

Planning policies are calculated with and without considering map uncertainties. The map is selected such that the mean values allow to calculate a feasible solution to the destination. As a consequence, the deterministic planning policy will correspond to the policy that is obtained using the Floyd-Warshall-based path planning approach from [START_REF] Kahveci | Optimal static soaring of UAVs using vehicle routing with time windows[END_REF]. This generates a comparison to the state of the art.

The simulations are performed using the simplified glider model from Section 3.2 with the polar curve parameters of the Cularis glider, as detailed in Appendix A. A one-shot planning is performed before the simulation and the map does not change while the vehicle is flying.

Therefore, the validity of a calculated policy holds for the entire flight and it suffices to calculate the policy once.

Updraft Map and Meteorological Parameters

Both the mean updraft parameters and the uncertainties are listed in Table 3 For each of the two maps, 100 true world parameters are sampled from the normal distribution defined by the respective map parameters. Note that the true world parameters are sampled such that there exists at least one feasible solution to the target. Then, the simulation is run with the probabilistic planning solution and the deterministic planning solution one time per sample. This enables to collect statistical information w.r.t. the defined maps.

Path Planner Settings:

The parameter k 1 is set to 0.01. The penalty for an outlanding, i.e.

the parameter k 0 is selected to be 5000 s to significantly penalize outlandings in the decisionmaking process. The granularities for the state discretization parameters m and n are 5 and 3 respectively, and the confidence interval for the ellipses is chosen at 99 %.

Meteorological parameters are required to solve for the planning horizon N , as previously defined in Equation (3.47). More specifically, both the average spacing between updrafts and the climb rates are required for calculating the average time a transit to thermal action requires. Regarding the spacing between updrafts, different observations can be found in the literature. While both the authors of [START_REF] Lindemann | Some characteristics of thermal convection as measured by a powered sailplane[END_REF] and [START_REF] Konovalov | Thermals in the sub-cloud layer of the atmosphere[END_REF] measured updraft spacings of 2.5 km, the author of [START_REF] Johnson | Measurements of airplane sink rates between thermals[END_REF] observed updraft spacings of 10 km. For the calculation of the planning horizon, we refer to the smaller value which will also hold for the higher updraft spacing values of 10 km. Typical updraft climb rates for cumulus soaring conditions are around w = 2.5 m/s [START_REF] Müller | Meteorologische Aspekte des Streckensegelfluges: Thermik von A bis Z[END_REF]. For the Cularis UAV, the MacCready airspeed corresponding to a climb rate of 2.5 m/s, i.e. V M C (w = 2.5 m/s) is 16.1 m/s at a sink rate of v si =1.53 m/s. Therefore, one can calculate a representative duration for an interthermal glide and climb sequence t aα

t aα = t g + t c = 2500 16.1 [s] tg=155 [s] + 1.53 t g 2.5 [s] tc=95 [s] = 250 [s].
Recall that updrafts last around 20 min [START_REF] Allen | Updraft model for development of autonomous soaring uninhabited air vehicles[END_REF]. This defines the maximum validity of the map t max = 1200 s, and one can determine the planning horizon

N = 1200 250 + 1 = 6,
where • denotes the operator for rounding up to the next integer value.

Evaluation: For evaluating the performance of both the deterministic and the probabilistic planning policies, both the mean realized Cost-To-Go (CTG) to the target, i.e. C (s 1 ) and the corresponding standard deviations are calculated for the two uncertainty levels (i.e. for the two maps). Recall that the CTG is defined as the sum of the overall flight time and the additive cost term Ψ as introduced in Equation (3.46), in the case of an outlanding. The statistics in terms of both mean values and standard deviations over all 100 MC runs are listed in Table 3 The mean cost of the probabilistic approach is approximately 1.5 times lower than its deterministic counterpart for the first map, i.e. at low uncertainties. At the increased uncertainty level, the probabilistic planning policy outperforms the deterministic approach by factor 4. The cost is optimized by selecting the seemingly fastest way. This is the trajectory from the origin over updraft #1 to the destination. Updraft #1 lies in proximity to the direct ground course.

Thus, only a tiny detour has to be tolerated which minimizes the flight time when solely considering the mean of the estimated parameters. However, if the vanish time of updraft #1 does not allow to climb back to an altitude from which the destination can be reached without undercutting the minimum altitude z min , the vehicle has no other option than to land out since the remaining updrafts do not lie in the vehicle's gliding range or have already vanished. Also note that the deterministic planner only provides one planning solution. That is, no back-up plans are calculated since the map is assumed to be perfectly known. In comparison to the deterministic plan, the probabilistic plan suggests a more conservative trajectory, initially guiding the vehicle via cloud #2. This plan increases the additional flight time since approaching the updraft #2 involves a detour, compared to the approach over updraft #1. However, an overestimation of the related vanish time still leaves the chance open to regain altitude by using updraft #3 in case the destination is not reachable after terminating the climb in updraft #2.

While the deterministic planner does not provide back-up plans, the probabilistic planning solution incorporates the necessitated back-ups for adapting the flight plan according to the true action outcomes.

A similar result is observed in Figure 3.19 where the flight trajectories are illustrated for the high map uncertainty case. Note that the trajectories of Figure 3.18a are identical to those in 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 In comparison, the probabilistic planner suggests a plan which is adapted to the higher level of map uncertainties and initially guides the vehicle over updraft #4, still leaving the option of regaining altitude using updraft #2 and updraft #3 respectively.

•10 4 -2,000 0 2,000 #1 #2
#3 #4 y [m] x [m] Updrafts Start Destination Outlanding site (a) Deterministic 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 •10 4 -2,000 0 2,000 #1 #2 #3
A further metric to compare the two approaches is the mean speed to the destination VD for successful missions, as illustrated in Figure 3.20. This parameter reflects the trade-off or cost that has to be tolerated in order to achieve less outlandings. As expected, the deterministic planning solution results in slightly higher ground speeds (approximately 1.1% and 2.5%). Note that the flight time increase induced by the probabilistic planner at high map uncertainties is due to the course deviation, when selecting updraft #4 as first action (see Fig- 3.2 and Figure 3.16 respectively, it is purposely plotted here to demonstrate that the suggested path planner works as expected, i.e. outlanding sites of the probabilistic planner generally lie in the proximity of the destination compared to the deterministic planner.

More specifically, the average remaining distance of the probabilistic path planner is approximately 80% less than for the deterministic path planner at both uncertainty levels.

Computational Cost

The time complexity of the path planning algorithm consists of both the complexity of the tree construction algorithm and the complexity of the backwards induction algorithm. Thus, the time complexity Γ of the path planning is defined by

Γ = O N |S| 2 |A| Γ bi + O |A||m × n| N -1 Γc , (3.48)
where the number of states is

|S| = N -1 k=1 |S k-1 ||m × n||A| k-1 and |S 0 | = 1. Note that Equa- tion (3.48
) is a weak definition since the time complexity generally refers to primitive computing operations. While the "min" function, as called by the backwards induction (Algorithm 1) can be considered as a primitive operation, the "calculate cost and probabilities" function which is defined as the primitive function of the tree construction (Algorithm 2) requires more elementary operations, when solving the integrals for both the probability and the expected value calculation.

Another crucial property is the runtime. Though the runtime depends on the code efficiency and the machine on which the algorithm is executed, it is interesting to show typical values, especially when the target hardware which is designated to be carried on-board the UAV features a similar processor as the computer on which the simulations are conducted. For the tree construction, the runtime depends on the number of actions |A|, the discretization settings m and n, as well as the planning horizon N . More specifically, the runtime depends on the number of function calls for "calculate cost and probabilities", i.e. on the number of generated states. However, it is challenging to generalize the number of actions accessible from each node and thus to analytically quantify the average or maximum number of states which translate in an average or worst-case runtime. This is because the number of actions that can be executed from a state will generally decrease with the layer which sensitively impacts the cumulated number of states and so the runtime. To understand this action reduction, consider the tree construction starting from layer one onwards. At this particular time instance, all updrafts in the map are still active and a majority of updrafts can be assumed accessible.

However, in the succeeding layers, time has elapsed and therefore the remaining lifespan of the updrafts has decreased which in turn reduces the number of applicable actions and thus the number of states.

In this work, we rely on the Monte Carlo method to identify both the average and the maximum (worst-case) number of states. This is done by generating a number of realistic inputs to the algorithm for capturing the statistics.

The number of states in the tree depends on the updraft density, the map size, both vanish times and updraft strenghts, and the state discretization parameters m and n respectively. Due to the circling trajectory of the UAV glider in thermal centering mode, the updraft map is best represented by a disk with a radius r d around the vehicle's current position where r d stands for the detection range of the vision sensor. Recall from Section 3.5.5, that course deviations up to ±90 • w.r.t. the target are tolerated when approaching updrafts. Therefore, the map features an area of 0.5πr 2 d (see Figure 3.22). Supposing both a perfect image processor and updraft estimator, the map contains all updrafts within the range r < r d which is a worst-case scenario for the computational load. One might argue that the map can contain more updrafts due to mismatched data association. However, the updraft position estimation is separated from the updraft strength and vanish time estimation, as was explained in Section 2.4.2. Further, the map only includes updrafts for which the vanish time and updraft strength estimation has been triggered. The trigger requires a consecutive series of matched vision measurements such that the threshold covariance γ is undercut. This procedure inherently reduces the likelihood of false updrafts in the map.

The detection range of the image processing sensor r d is selected to be 5000 m. The updraft density ρ u can be defined as function of the updraft spacing d u . For a regular distribution, we can rely on the density formula of a hexagonal grid. Thereby, the updraft density is defined

ρ u = 4 √ 3d 2 u . (3.49)
By means of the updraft density, we can calculate the number of updrafts that lie in a given area of size A

n u = ρ u A. (3.50)
This yields n u = 0.5ρ u πr 2 d . updrafts for the half disk. The parameters used to generate sample trees are listed in Table 3.3.

Parameter

Value Unit is likewise shown. Note that this worst case holds under the condition that all updrafts are accessible from every state. As expected, both the maximum and average number of states at each layer, resulting from the Monte Carlo study lie significantly lower than for the theoretical worst-case. Also note that for increasing decision epochs, the number of states decreases. This is due to the fact that only a fraction of the updrafts last long enough to enable "long" routes, i.e. routes where the vehicle visits more than 4 updrafts before transitioning to the target. To illustrate this effect independently of the discretization parameters m and n, the ratio between the number of states obtained from the MC study and the theoretical worst-case is illustrated in Figure 3.24. The reduction of states becomes especially significant in the third and fourth layer, where the Monte Carlo study predicts only 25% and 3% of the theoretical worst-case respectively. The corresponding runtimes for both the tree construction and the value iteration are listed in Table 3.4. Note that the algorithms were programmed in Matlab®, mex-compiled, and executed on an Intel Core i7 processor. Also note that miniaturized COTS computing systems with equivalent specifications and dimensions compatible with small UAV gliders are readily on the market. As expected, the runtime required for the tree construction is significantly higher than the runtime for the backwards induction. More specifically, it requires about 16 times more to calculate the tree than to solve the MDP. However, for realistic input settings (number of updrafts) and a worst-case scenario where all updrafts are identified by the estimator, an average runtime of approximately 6 s is still considered to be fast enough regarding the fact that the relative motion between the vehicle and the observed clouds is slow -resulting in slow changes of the updraft maps which does not require for high re-planning rates. Additionally, the code was not optimized w.r.t. runtime which means the runtime can still be reduced, for example through vectorization. Additionally, a study is conducted, demonstrating that the proposed probabilistic planner is shown to be implementable on a small computer, enabling an on-board execution using Commercials Off-The-Shelf (COTS) electronic devices.

d u 2.5 [km] r d 5 [km] n u 15 [-] t v U(0, 20) [min] w N (2.5, 0.5) [m/s] m 5 [-] n 3 [-] N 6 [-]

System Architectures

In order to highlight the benefits of both the remote updraft estimation and the probabilistic path planner, the simulations in this chapter are repeated for different system architectures. That is, we define three systems that consist of different constellations of estimation and planning components. The first one aims to reproduce the state-of-the-art performance in autonomous cross-country soaring. Therefore, the glider flies without any remote updraft sensing which reproduces flight performances that are obtained using the techniques of [3],

i.e. only in-updraft estimation. The two other architectures rely on the vision-based remote updraft estimation from Chapter 2 where the second system architecture is based on the deterministic path planner and the third one on the probabilistic path planner. This allows for both a rigorous and sequential evaluation of the elaborated navigation and planning methods.

Thermal Centering Control For all three system architectures, in-updraft identification is realized by monitoring the vehicle's energy state. Note that the energy state is defined by both the first and the second derivative of the total energy (E = mgh + 1/2mV 2 a ), i.e. both total energy rate Ė and acceleration Ë. A mode logic, suggested by the authors of [2], outputs an "in-updraft flag" as soon as the vehicle's energy state passes a predefined threshold.

Further, all three architectures rely on the same thermal centering control law which is inspired by the work presented in [START_REF] Reichmann | Cross-Country Soaring[END_REF]. The centering control law builds on three rules of thumb in order to generate bank angle controls that are function of the vehicle's energy state 1. As climb improves, decrease the bank angle, 2. As climb deteriorates, increase the bank angle, 3. If climb remains constant, maintain a constant bank angle.

These three rules are implemented to generate a commanded heading derivative

Ψcmd = Ψss + K tc Ë, (4.1)
where Ψss is a steady state feed forward contribution which is function of the airspeed and some selected steady state radius r ss such that

Ψss = V a /r ss (4.2)
and K tc is the control law gain. Note that the turn rate Equation (4.2) is an approximation, holding for small flight path angles which is mostly the case for glider planes. The bank angle command for thermal centering is then defined assuming a coordinated turn by

Φ cmd = arctan V a g Ψcmd . (4.3)
System Architectures The three system architectures are detailed as follows:

1. Blind: The first architecture is non-vision based, i.e "blind". Updrafts are thus not remotely sensed and the vehicle "blindly" flies towards the target, switching to thermal centering mode whenever the presence of an updraft is identified, given the vehicle's energy state. For this architecture, airspeed is always set to the value of the best glide

Monte Carlo Simulations

For studying the performances of the three system architectures, a series of Monte Carlo simulations is conducted.

Mission

The benchmark mission of the Monte Carlo Simulation (MCS) is a 25 km cross-country soaring flight where updrafts are randomly located within a 25 km× 7 km weather tile around the direct line between the start point and the target (Figure 4.4). The vehicle is released at the origin at an altitude of 1100 m in thermal centering mode and its task is to fly to the destination. Recall that the vehicle is not allowed to undercut a minimum altitude, selected to be 500 m. Thus, mission success is whenever the vehicle reaches the destination. The simulation is stopped as soon as an outlanding occurs, i.e. the vehicle's altitude undercuts the minimum altitude or whenever the vehicle reaches the target.

Monte Carlo Setups

In cross-country soaring, the spacing between updrafts has a major impact on the mission success. Since this parameter can significantly vary, we define two different levels of updraft densities to increase the validity of the simulation results and also to investigate the robustness of the individual architectures. A further crucial parameter for the two vision-based system architectures is the image processing noise, as introduced in Equation (2.25). Therefore, as for the updraft density, two levels of image processing noise are defined. This generates four MC setups, i.e. the permutations of low/high updraft density and low/high image processing noise. For each of these four setups, a total of 100 MC simulation runs is conducted for both the second and the third system architecture. Since the performance of the "blind" system architecture is independent from the image processing noise, MC simulations for this architecture are solely conducted at the low and the high updraft density settings. Note that the scenarios are feasible, i.e. there exists at least one solution to the target.

Both the updraft spacings d u and the image processing noise levels Σ y LR defining the four MC setups are listed in Table 4.1. The values for the updraft spacings are taken from the literature, where the lower spacing refers to the findings of [START_REF] Lindemann | Some characteristics of thermal convection as measured by a powered sailplane[END_REF][START_REF] Konovalov | Thermals in the sub-cloud layer of the atmosphere[END_REF], and the upper updraft spacing to ones of [START_REF] Johnson | Measurements of airplane sink rates between thermals[END_REF]. Regarding the setting of the image processing noise, the default level from Chapter 2

is selected as low level, and twice its magnitude is selected to define deteriorated conditions.

Parameter

n u = ρ u A T = 4 √ 3d 2 u A T , (4.4

Outlanding Statistics

Since the main purpose of the vision-based soaring approach is an augmentation of the reliability of autonomous cross-country soaring systems, the cumulated number of outlandings constitutes the most important figure of merit. In the same figure, the number of outlandings at high image processing noise are listed. Since the blind system architecture is insensitive to image processing noise, it results again in 86 outlandings over all runs. Evidently, the increased image processing noise causes an increased number of outlandings for the two vision-based system architectures due to a degradation of the remote updraft estimation. More precisely, the number of outlandings for the deterministic planner increases by 3 while the outlandings caused by the probabilistic planner increase by 7. As for the low measurement noise, the trend remains, i.e. the probabilistic architecture causes again less outlandings, compared to the deterministic planner.

Figure 4.5b depicts the number of outlandings for low updraft density at both image processing noise levels. Interestingly, the "blind" system architecture not only degrades, but features a mission fail rate of 99%. At low updraft density and low image processing noise, the visionbased approaches have success rates of 54% for the deterministic and 58% for the probabilistic planner, i.e. the probabilistic planner still outperforms the deterministic one. For low updraft density and high measurement noise the two vision based approaches feature comparable performances. That is, the probabilistic planner yields one outlanding less than its deterministic counterpart. Figure 4.5 further reflects how the updraft spacing/density impacts the success rate. While the "blind" performance features poor success rates at both density levels, and as such is relatively insensitive to the regarded updraft density, the vision based approaches remain sensitive to both the updraft density and the image processing noise with a success rate variatin from 47% to 92% for the deterministic, and from 48% to 100% for the probabilistic path planner respectively.

Analysis and Comparison of Flight Trajectories

In this section, a representative MC simulation sample is studied to illustrate the different functionalities of the three system architectures. Using the same meteorological setup with an updraft spacing of 2.5 km, and the same noise seed, the flight is repeated three times, each time with another system architecture. The trajectories of the probabilistic and the deterministic architecture coincide up to the split point. Therefore, both planners have exactly the same state estimates at their disposal until this point. However, both the estimates of the updraft strength and vanish time of the fourth updraft (see Figures 4.8e to 4.8h) are of poor quality. More specifically, the vanish time is overestimated by approximately 400 s. This overestimation and the target-proximity of the updraft explain the choice of the deterministic planner. While approaching the updraft, the estimate converges until the point where the path planner judges that an approach is not beneficial. This instance is illustrated by a gray point (stopp approach point) in Figures 4.6 to land out since a potential approach to any other updraft would violate the minimum altitude constraint when flying with the MacCready speed. Note that this outlanding could have been avoided by a reduction of the MacCready airspeed. Thereby, this example also reveals a potential refinement of the proposed path planner that can be obtained by optimizing both airspeed and updraft selection in parallel. Contrary to the deterministic planner, the probabilistic one exploits the uncertainty information of the updraft estimates and continues to guide the glider towards updraft #3 which bears only a neglectable risk of a missjudged approach. Figures 4.8g

Flight Time Statistics

As a second figure of merit, the flight time for successful flights is studied. The flight time statistics for all four MC setups are illustrated and denoted in Figure 4.9 and Table 4 and 4.4b. Note that the blind system architecture is insensitive to image processing noise and thus features equal performances at both noise levels. Independently from the image procesing noise, at high updraft density, the deterministic architecture yields the lowest mean flight time to the target with 49.42 and 50.81 min respectively. One could expect the blind system architecture to yield the lowest mean flight time because no detours are flown. This is however not the case since the vehicle does not fly with the MacCready airspeed but with the airspeed of the best glide ratio which is considerably lower. At both image processing noise levels, the probabilistic architecture results in a mean flight time increase of approximately 1 min, compared to the deterministic architecture thereby reflecting the more conservative planning approach.

The flight time increase of around 1.5% can however be tolerated regarding the significantly higher success rate of the probabilistic architecture. Notice that at increased measurement noise, the mean flight time of both vision-based approaches increases. This is because the degradation of the estimation causes more missed approaches to updrafts. That is, if the vanish time of an updraft is highly overestimated, the corresponding flight leg adds to the flight time but the vehicle does not gain energy. The standard deviations of the flight times for both vision-based architectures are slighty higher than for the blind system architecture. This effect is mainly introduced by the randomness of the updraft positions that can create significant flight detours. levels, the vision-based approaches result in comparable mean flight times. This is because the decision-making process is less complex than for the high updraft density case, increasing the likelihood of both planners to make similar decisions. Interestingly enough, the flight times at low updraft density are lower than at high updraft density. This is explained by the fact that less updrafts are visited. Therefore, less time is spent in thermal centering mode which evidently reduces the flight time.

Analysis of the Mean Altitude

A further metric that indicates the differences between the three system architectures is the mean altitude at each time step over all simulation runs, i.e for both failed and successful flights, as illustrated in Figure 4.10 for all three architectures. Note that the peak at approx. 150 s and 1200 m is the point where the glider leaves the initial updraft at the origin of the coordinate system. First, the figure illustrates that the blind system architecture features a significant variation in altitude towards the target. The mean altitude is dominated by failed missions, causing outlandings after approximately 1500 s. This phenomenon could already be observed in the flight trajectory study, as illustrated on Figure 4.7a. In comparison, the two vision-based architectures maintain a higher average altitude. Furthermore, the altitude of the vision-based architectures remains almost constant between the first leg and the last leg. Also note that the more conservative planning of the probabilistic planner is reflected in a slightly higher mean altitude, compared to the deterministic planner.

Cost Function Variations

In addition to the MC simulation and sensitivity study, further simulations were conducted to illustrate the impact of varying cost function settings of the proposed path planning method.

Recall from Chapter 3 that the cost for transition to target actions is function of the two parameters k 0 and k 1 . The first parameter serves to penalize non-feasible paths in general, and The impact of the parameter setting can also be observed by counting the number of updrafts that the vehicle approaches. In the first two cases, the vehicle approaches 5 updrafts (see the small dots that lie on the trajectories). Conversely, in the case of the most aggressive setting, only 2 updrafts are approached. Note that the flight was executed using system architecture 3, i.e. using vision and the probabilistic planner, but the same effects can be observed using the deterministic planner.

Summary

In this chapter, an overall evaluation of the proposed navigation and path planning methods for autonomous cross-country soaring is presented. For enabling both a detailed evaluation of the presented methods and highlighting the benefits compared to state-of-the-art approaches, Chapter 5

Conclusion

In this chapter, the thesis' contributions are summarized, and directions for future research in this field are discussed.

Autonomous soaring for fixed-wing UAVs is a relatively new field of research, and a technique that offers significant potential to increase both the flight time and the range of unmanned aircraft. However, a key challenge in autonomous soaring lies in the identification of updrafts and in the generation of efficient flight paths.

Therefore, this thesis has proposed a framework for Vision-Based Autonomous Cross-Country Soaring (VBACCS). The originality of this work lies in the use of vision for remotely estimating the position and parameters of sub-cumulus thermal updrafts, thereby mimicking the behavior of human glider pilots. While the realization of a VBACCS system necessitates progress in various fields, the focus was put on the remote estimation and path planning part. Although the primary purpose of this thesis is autonomous soaring, some of the contributions find application beyond that field, e.g. in meteorological research.

Thesis Contributions

Vision-based Framework for Autonomous Cross-Country Soaring

The first contribution of this thesis is a new framework for reliable autonomous cross-country soaring flight of a small glider UAV that relies on vision for remotely estimating updrafts that are capped by cumulus clouds. The idea for the framework is inspired by the behavior of human glider pilots who scan the environment for cumulus clouds and use this information in their decision-making. Parts of this concept are realized by proposing algorithms for visionbased updraft estimation and map-based path planning complemented by an evaluation in a simulation environment.

Vision-based Remote Updraft Estimation

In this thesis, an estimator for remotely estimating updrafts is developed. Regarding the state of the art before this thesis, the estimator provides an original solution to the remote updraft estimation problem.

The UKF is identified as suitable candidate due to its capability of providing consistent state estimates for non-linear systems. The proposed estimator features a hybrid architecture, in which the estimation of updraft positions and updraft parameters are performed separately in order to cope with varying initialization dynamics.

Both vision measurements and the vehicle's attitude and position are taken as filter inputs in order to calculate the position of cumulus clouds in the first filter. The phenomenon of a uniform cloud base altitude was exploited by estimating cloud positions in one filter, rather than in individual filters in order to increase estimation convergence. Furthermore, a virtual cloud width sensor was designed. The sensor transforms vision measurements and the cloud's estimated position from the first filter into virtual measurements of the current cloud width. The technique of Unscented Transform (UT) is applied to reliably propagate measurement uncertainties through the system equations of the sensor. The proposed cloud width sensor can also be used for measurement campaigns in meteorological research.

In the second filter, the measurements of cloud widths are used for tracking their evolution over time. True measurements of cloud width evolutions were used for identifying a simple function that describes the cloud size dynamics. This function is relied on to predict both the vanish time and the updraft strength related to the cumulus cloud. MCSs were conducted in order to study both the convergence and consistency characteristics of the proposed estimator using a realistic observer soaring flight trajectory. The results show the capacity of the proposed filter to provide both convergent and consistent estimates in spite of adverse observation conditions.

Robust Path Planning For Map-Based Autonomous Cross-Country Soaring

Map-based path planning for autonomous cross-country soaring deals with the calculation of a series of updrafts that the glider has to approach for regaining altitude before it can reach its destination. However, under real world conditions, the updraft maps are imperfect due to the presence of measurement errors. Therefore, a robust path planner for map-based autonomous cross-country soaring was developed. The originality of the path planner is that it explicitly handles the trade-off between reducing flight times and a reduction of the outlanding risk by considering updraft map uncertainties.

The framework of MDPs was identified as suitable candidate. The path planning problem was formalized in a computationally efficient way by using a simple flight dynamics model. Further, a final glide guidance law was desinged to minimize the flight time of the final leg, while respecting constraints on the minimum altitude. Using the flight dynamics model and the technique of Unscented Transform (UT), the cross-country soaring problem can be modeled as a finite horizon MDP. A cost function based on a glider-competition rule was suggested for the calculation of transition costs and backwards induction was applied to determine the seeked optimal policy. A series of simulations illustrate the benefits of the proposed path planner in comparison to a state-of-the-art deterministic planner. Additionally, a study of the execution time of a prototype implementation is provided, indicating its implementability on-board a small scale UAV glider, equipped with low-cost computer devices.

Evaluation of the Overall Vision-Based Autonomous Cross-Country Soaring (VBACCS) System

To assess the joint performance of both the vision-based remote updraft estimation and the robust path planner, a series of Monte Carlo Simulation (MCS) was conducted using a realistic 6 DOF flight simulator including low-level flight controllers and a thermal centering controller.

Three different system architectures for autonomous cross-country soaring were defined and compared. One of these reproduces the state of the art in autonomous cross-country soaring, prior to this thesis. Thereby, performance of the proposed navigation and planning algorithms is not only individually assessed, but also compared. Representative meteorological parameter settings and varying image processing noise levels were used in the simulations in order to also study the sensibility of the proposed system with respect to crucial input parameters. The simulation results confirmed the significant potential of the proposed VBACCS system to contribute to a reduction of outlandings, while yielding only minor increases in flight time.

Limitations and Suggested Future Work

Working over a long time on a complex subject generates not only solutions, but also reveals questions and even more challenges. This section discusses the limitations of the presented work and suggests directions for future work.

• Horizontal Wind: In the presence of horizontal wind, the position of clouds and updrafts might not be stationary but change with time. To consider wind in the estimation process, the glider could be equipped with a wind estimator and the updraft estimator could rely on the wind estimates in the state prediction function. Another possibility would be to estimate both wind and cloud positions in one filter. Therefore, airdata (e.g. airspeed and angle of attack) and vision measurements would have to be fused. Moreover, models should be integrated that incorporate possible updraft tilting and bending due to horizontal wind. In the same vein, the effects of horizontal wind will have to be integrated in the state propagation model of the suggested path planner.

• Attitude Uncertainties: In this thesis, both attitude and position estimates of the vehicle were assumed to be of high quality. Degraded attitude estimates however propagate into the updraft estimation and can sensitively impact the estimation performance. An extension of the proposed estimator would therefore be to incorporate the uncertainty of attitude and vehicle position estimates in the corresponding system equations.

• Position Uncertainties: In this work, updraft position uncertainties were not considered in the path planning. However, especially in the presence of horizontal wind, increased updraft position uncertainties can be expected. This can impact the seek-time under an approached cloud. Further research should therefore concentrate on how these position uncertainties can efficiently be integrated into the proposed path planner, while keeping the problem tractable.

• Image-In-The-Loop-Simulation: Another important direction of research will include the evaluation of the estimator with images in the loop. In a first step, synthetic images that reproduce realistic lightning and occlusion effects could be used. This requires a suited image generator. A possible direction would be to rely on images from a computer game flight simulator with a high degree of realism. Both the estimator and the path planner could thereby be tested in "closed-loop".

Going one step further, flight data recordings and vision measurements should be collected for evaluating the estimator with true inputs. These inputs can for an example be collected using a manned motor glider plane.

• Joint Airspeed and Updraft Selection: In the proposed path planner, updraft and airspeed selection were separated. This procedure allows to keep the planning problem tractable, but also has its shortcomings, as was shown in Chapter 4. Therefore, future research should concentrate on how to efficiently combine the optimization of both the updraft and the airspeed selection.

• A Priori Knowledge And Fusion Of Soaring Techniques: Complementary to information on remotely located updrafts, it is also beneficial to incorporate geographical a priori knowledge in the path planning task. It can for an example be unfavorable to fly over open waters such as lakes or seas, since it is more likely to encounter temperature gradients and thus updrafts over land. In the same vein, it could be promising to jointly benefit from different soaring techniques (e.g. ridge soaring, dynamic soaring etc.) during one autonomous soaring flight. Therefore, the path planner could be extended to rely on geographical information together with the perceived environment in order to decide when to exploit which energy source and to set the corresponding mode of operation.

• Adaptive Cost Function: The proposed path planner relies on a cost function that incorporates the notion of an outlanding. Throughout this thesis, the cost function parameters for outlandings were kept constant during flight. Meteorological conditions in terms of both the updraft density and strength can vary during the course of the day. While updrafts are comparably weak in the morning and the early evening, stronger updrafts are encountered during the afternoon. For long distance flights, lasting several hours, the cost function of the proposed path planner could therefore be adapted based upon the current meteorological conditions. To give an example, if it is likely that updrafts will occur on the future flight path, one can include this information by changing the cost function to be less conservative.

• Target Hardware Implementation and Flight Testing: Another important aspect is the implementation of the proposed algorithms on the designated hardware components.

While prototypes of the systems were programmed for simulation purposes, code efficiency plays an important role for embedded systems, especially when relying on the computationally demanding framework of MDPs. Therefore, these prototypes should be optimized. Once the navigation and path planning components have been implemented, the entire system can be evaluated in real flight tests. Flight tests require both suitable meteorological conditions, as well as flight clearances, since the UAV glider will most likely fly outside the line of sight of the operator. Therefore, it could be interesting to conduct these flight test campaigns over uninhabited and desert-like regions such as the Mojave desert in California or the Namibia desert, thereby facilitating flight clearances while also benefiting from advantageous meteorological conditions that help to evaluate the algorithms.

A.1.1. Aerodynamic Forces and Moments

The aerodynamic forces F a and moments M a (in the body frame) that act on the vehicle are defined as Note that all angles are given in rad.

F a = R b←s qS      C D C Y C L      , M a =
The airspeed controller Equation (B.2) generates appropriate short period mode damping around values of 0.7. Both the airspeed error V a,cmd -V a and the integrated airspeed error as well as the angle of attack α, the pitch angle Θ, and the pitch rate q are fed back.

δ e = K Va
0.013 0.006 0.002 -0.219 0.005 The conditional expectation of the transformed variable x is defined by The parameter c 2 is obtained in a similar manner. Integration over the inner variable x1 results in Rearranging the terms, we obtain Calculating the integrals of c 1 and c 2 for every subregion, we obtain the conditional expectation in the transformed coordinate system which is then finally transformed into the original frame.

            V a,cmd -V a V a,cmd -V a dt α trim -α Θ trim -Θ q             qref q (B.
E x|x ∈ B = 1 P x ∈ B 3 i=1            1 2π|Λ| 1
c 2 = 1 2 √ π √ 2 √ λ 1 x2,u x2,l x2
c 2 = 1 √ 2 √ π √ λ 2 e x2 2λ 2 erf α u x2 + β u √ 2λ 1 -erf α l x2 + β l √ 2λ 1 x 2,u x 2,l + √ λ 2 2π √ λ 1   e - β 2

  Reference scenario . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 2.15 Normalised Estimation Error Squared (NEES) and Root Mean Square Error (RMSE) for the cloud position tracking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 2.16 NEES for updraft vanish time and updraft strength . . . . . . . . . . . . . . . . 48 2.17 RMSE for updraft vanish time and updraft strength . . . . . . . . . . . . . . . . 48 2.18 Estimation error depending on the number of tracked clouds . . . . . . . . . . . 50 3.1 Path planning problem for map-based cross-country soaring . . . . . . . . . . . 54 3.2 Interthermal glide . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.3 Glider kinematics in the vertical plane during glide . . . . . . . . . . . . . . . . . 56 3.4 MacCready airspeeds as function of the climb rate, calculated for the Cularis UAV 57 3.5 V a as function of the excess altitude p z -z min and the horizontal distance q D to the destination . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 3.6 The effect of p z,e on outlandings . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 3.7 Probability density of the uncertain terminal state . . . . . . . . . . . . . . . . . 65 3.8 Confidence ellipse of s ep with related bounding box . . . . . . . . . . . . . . . . 72 3.9 Discretized saturated bounding box . . . . . . . . . . . . . . . . . . . . . . . . . 73 3.10 Bin-related probabilities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 3.11 Discrete exit states with related bin probabilities . . . . . . . . . . . . . . . . . . 77 3.12 Side view of the discretized exit states . . . . . . . . . . . . . . . . . . . . . . . . 77 3.13 Cost calculation for transit-to-target actions . . . . . . . . . . . . . . . . . . . . . . 78 3.14 Illustration of the decision tree . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 3.15 Reference scenario, xy-plane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 3.16 Cost-To-Go (CTG) over sample . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 3.17 Cumulated number of outlandings . . . . . . . . . . . . . . . . . . . . . . . . . . 86 3.18 2 D-Flight trajectories and outlanding sites at low map uncertainties . . . . . . . 87 3.19 2 D-Flight trajectories and outlanding sites at high map uncertainties . . . . . . 88 3.20 Average speed to the destination . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 3.21 Average remaining distance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 3.22 Half-disk-shaped updraft map . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 3.23 Nodes per layer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 3.24 State number ratio between the theoretical worst-case scenario and the Monte Carlo results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 4.1 Cularis electric glider (image taken from http://multiplex-rc.de/) . . . . . . 96

AlphanumericbLδ 2 Σ

 2 , b R Left/Right bearing to the cloud feature points, m b Wingspan, m w Climb rate of the glider, m/s d Cloud size evolution vector, m, m/s, m/s 2 d u Updraft spacing, km f Focal length, m p Vehicle position vector (Inertial frame), m p c Cloud position vector, m p u Updraft position vector (inertial frame), m q d Horizontal distance between the vehicle and the destination position, m q r Horizontal distance between the outlanding site and the destination position, m S x , S y Conversion factors from meter to pixel, px/m Ground speed, m/s V M C MacCready airspeed, m/s v si Sink rate, m/s w p Core updraft velocity, m/s w p,z Core updraft velocity averaged over altitude, m/s w Convective velocity scale, m/s w u Updraft velocity, m/s wu,rz Average updraft velocity of all updrafts with respect to both the updraft radius and the altitude, m/s z Altitude, m z i Convective mixing layer thickness, m z cb Cloud base altitude, m Greek letters Angle between the heading to the cloud base center and the feature points of the cloud, rad γ Flight path angle, rad κ Coefficient relating the updraft velocity to the cloud diameter amplitude, 1/s ν Zero mean Gaussian measurement noise Ψ Heading, rad Ψ Penalty for outlandings, s ρ u Updraft density, 1/km Diagonal matrix with standard deviations τ c Lifespan of a cloud, s Chapter 1

Figure 1 .

 1 Figure 1.1 illustrates the flight trajectory of an Albatross during dynamic soaring. Nature has motivated researchers in bringing this capacity to UAVs [1, 8]. Although dynamic soaring for

Figure 1 . 2 :

 12 Figure 1.2: Illustration of ridge soaring (http://www.aerospace.org)

Figure 1 .

 1 4 illustrates the glider plane used by Edwards. The UAV could perform a fully autonomous cross-country soaring flight over a distance of approximately 50 km. The corresponding flight trajectory is depicted in Figure 1.5.

Figure 1 . 6 :

 16 Figure 1.6: Cumulus soaring conditions in Pokweni, Namibia (http://www.gliding.ch)

Figure 1 . 7 :

 17 Figure 1.7: Flight log of a 1400 km FAI triangle over Namibia (data and figure taken from http://www.onlinecontest.org)

  Figure 1.9: A vision-based framework for reliable autonomous cross-country soaring

Figure 2 . 1 :

 21 Figure 2.1: Setup between glider, camera and cloud

. 12 )

 12 Cloud position: With the definition of the updraft position from Equation (2.8), both cloud and updraft position coincide. Therefore, the following holds p c = p u and ṗc = ṗu .(2.13)

Figure 2 . 2 :

 22 Figure 2.2: Uniform cloud based altitude

Figure 2 . 3 :

 23 Figure 2.3: Field of cumulus humilis clouds

Figure 2 . 4 :Figure 2 . 5 :

 2425 Figure 2.4: Lifecycle of a cumulus

Figure 2 . 6 :

 26 Figure 2.6: Fitted and measured cloud width evolutions

Figure 2 . 7 :

 27 Figure 2.7: Pinhole camera model

  vision-based updraft estimation builds upon the fact that one can determine the updraft position, strength, and lifespan by observing cumulus clouds. As a minimum visual information, the two feature points y L and y R which are the endpoints of the cloud's base in the image are required. These points are shown in Figure 2.8.

Figure 2 . 8 :

 28 Figure 2.8: Feature points

Figure 2 . 9 :

 29 Figure 2.9: Image processing steps

Figure 2 . 10 :

 210 Figure 2.10: Geometry between y c and y LR

Figure 2 . 11 : 2 . 4 . 3 .

 211243 Figure 2.11: Structure of the estimator

Figure 2 . 12 :

 212 Figure 2.12: Spiral-shaped flight trajectory during thermal centering flight, [3]

  Vision measurements y c of the clouds center point are predicted from estimated cloud positions p c = [p T cxy z cb ] T using the camera model from Equation (2.22)

  Each time a new cloud is detected, both its estimated position and covariance have to be initialized. As the clouds lie on the horizontal plane which is defined by the cloud base altitude z cb , it is straightforward to calculate the initial cloud position as the plane-line intersection between the cloud base plane and the line defined by the bearing b v c

  Figure 2.14a. The vehicle centers an updraft at 950 m altitude and climbs to the cloud base at about 1200 m altitude. The flight trajectory is generated relying on the aircraft model and the flight control laws that are described in Appendices A and B. The vehicle circles around the origin of the xy frame, where the lowest circle in Figure 2.14a is slightly shifted with respect to the upper circles. This is because the vehicle is not yet centered around the updraft's core at the beginning of the simulation. During the climb, the updrafts of four clouds are estimated, as shown in Figure 2.14b. The rhombi in Figure 2.14 indicate the mean value of the true cloud positions. However, for each simulation run, the true cloud positions are varied, where the variation is specified by the mean and a given standard deviation.

  Figure 2.14: Reference scenario

Figure 2 .Figure 2 . 15 :

 2215 Figure 2.15 illustrates Monte Carlo results in terms of NEES and RMSE from 100 runs for the position estimation evaluation. The dashed lines in the upper graph (Figure 2.15a) represent the two sided acceptance interval for N = 100 runs and a state dimension of L = 9 for a probability concentration region of 95%. Note that the state dimension L is determined from 4 × 2 D cloud positions and the cloud base altitude z cb .The NEES remains within the acceptance bounds, i.e. the covariance of the position state is correctly estimated, despite the adverse observation conditions. The consistent position estimate provides an important result for the vanish time and updraft velicity estimation. This is because the measurement noise of the cloud diameter sensor from Equation (2.60) depends

Figures 2 .

 2 Figures 2.16 and 2.17 show the results of the second Monte Carlo simulation -dedicated to the evaluation of the updraft strength and vanish time estimation. The four solid lines in each of the graphs specify the results of the individual clouds. Note that the graphs illustrate the NEES and RMS of the vanish time and updraft velocity, obtained after propagating the estimated parabola parameters through the updraft cloud relation.

Figure 2 .Figure 2 . 16 :

 2216 Figure 2.16 illustrates the result in terms of consistency with the dashed lines defining the region of acceptance. For all clouds, the NEES mostly remains within the acceptance bounds, indicating that both the uncertainty of the updraft strength and vanish time are correctly estimated.Another interesting fact is that the initial RMSE, illustrated in Figure2.17b, decreases slowly, compared to the RMSE of the position estimation from Figure2.15b. This is due to the fact that

Figure 2 .Figure 2 . 18 :

 2218 Figure 2.18a illustrates the RMSE of the horizontal position p cxy of cloud #1. The updraft strength and vanish time estimates of that cloud are shown in Figures 2.18b and 2.18c. Note that in both graphs, the bold line represents the case where four clouds are tracked.As expected, the error of the horizontal position decreases faster as the position estimate of the cloud is updated even while no vision measurements corresponding to this cloud are available. This can be seen at t =25 s where the cloud disappears for the first time from the vision sensor. For the case where only one cloud is tracked, the horizontal position error remains constant until the glider has performed a full circle and the cloud reappears in the camera's field of view. This is in contrast to the case when observing four clouds. Even while the cloud number one is not in the field of view, its horizontal position error decreases. These varying estimation performances in terms of convergence time propagate into the time at which the vanish time and updraft strength is triggered. Therefore, consider Figures 2.18b and 2.18c. In the multi-cloud run, the updraft parameter estimation of cloud 1 is triggered at t =75 s which is approximately 50 s earlier than in the single-cloud case. In other terms, information on the updraft strength and vanish time is available earlier, and the path planning can compute trajectories based on more complete information of the vehicle's environment.

. 1 ,

 1 the cross-country soaring mission for the glider consists of flying to a destination p D in minimum time, using updrafts that arise during the flight. The destination is given by p D ∈ R 3 where [p Dx p Dy p Dz ] T . The altitude of the vehicle at the destination p Dz is a free parameter such that p Dz ≥ z min where z min is some minimum altitude. Therefore, the path planning can be seen as a constrained shortest path problem where the vehicle's vertical operational range is restricted to z min ≤ p z ≤ z cb . Recall from Figures 2.2 and 3.1 that z cb is the cloud base altitude. Undercutting the minimum altitude z min is considered as mission failure.

Figure 3 . 1 :

 31 Figure 3.1: Path planning problem for map-based cross-country soaring

Figure 3 Figure 3 . 2 :

 332 Figure 3.2.

Figure 3 .

 3 4 illustrates the relation between w and V M C for the polar curve of the glider from Appendix A.

Figure 3 . 4 :

 34 Figure 3.4: MacCready airspeeds as function of the climb rate, calculated for the Cularis UAV

Figure 3 . 5 :

 35 Figure 3.5: V a as function of the excess altitude p z -z min and the horizontal distance q D to the destination

Figure 3

 3 m 0 m 1200 m 400 s 3 m/s T ,

Figure 3 . 7 :

 37 Figure 3.7: Probability density of the uncertain terminal state

  k = N -1 and successively decrementing k by one, the value of all states is minimized according to V k (s) ← min a∈A s ∈S P(s |s, a) C(s |s, a) + V k+1 (s ) . (3.29)

  |s, a) C(s |s, a) + V k+1 (s ) .(3.30) Due to the reversed chronological optimization, this procedure is known as backward induction[START_REF] Puterman | Markov Decision Processes: Discrete Stochastic Dynamic Programming[END_REF][START_REF] Powell | Approximate Dynamic Programming: Solving the Curses of Dimensionality[END_REF]. The pseudo code for the backward induction procedure is noted in Algorithm 1. The asymptotic time complexity (noted O) of the backwards induction algorithm depends on both the number of the backups and the number of the states. The complexity of one backup is O |S| 2 |A| where | • | represents the number of elements contained in the respective set. Thus, the time complexity for a horizon-N backward induction is O N |S| 2 |A| .

Algorithm 1 : 1 .

 11 Backwards Inductionbackwardsinduction() Input : MDP with (S, A, P, C, K) Output: V (s) min. value π (s) optimal policy Initialize: V N (s N ) = 0 for k = N -1 to k = 0 do for each non-goal state s ∈ S k do V k (s) ← min a∈A s ∈S P(s |s, π(s)) C(s |s, π(s)) + V k+1 (s ) π k (s) ← arg min a∈A s ∈S P(s |s, π(s)) C(s |s, π(s)) + V k+1 (s ) MDP State Definition For the cross-country soaring problem, the set of states S, as defined in Section 3.4.2, contains all possible 4 D position-time states s, noted as s = (p x , p y , p z , t) . (3.31)

Figure 3 . 11 :Figure 3 . 12 :

 311312 Figure 3.11: Discrete exit states with related bin probabilities

Figure 3 . 13 :

 313 Figure 3.13: Cost calculation for transit-to-target actions

Algorithm 2 :

 2 Decision tree constructionsconstructdecisiontree() Input : MDP with (A,s 0 , m, n and horizon N ) Output: tree for k = 0 to k = N -1 do for each s ∈ S k do for each a ∈ A do calculate cost and probabilities of all m × n action outcomes

Figure 3 . 15 :

 315 Figure 3.15: Reference scenario, xy-plane

Fig- ure 3 .Figure 3 . 16 :Figure 3 . 17 :

 3316317 Figure 3.16: Cost-To-Go (CTG) over sample

Figure 3 .

 3 Figure 3.18: 2 D-Flight trajectories and outlanding sites at low map uncertainties

Figure 3 .

 3 Figure 3.19: 2 D-Flight trajectories and outlanding sites at high map uncertainties

ure 3 .Figure 3 . 20 :

 3320 Figure 3.20: Average speed to the destination

Figure 3 . 21 :

 321 Figure 3.21: Average remaining distance

Figure 3 . 22 :

 322 Figure 3.22: Half-disk-shaped updraft map

Figure 3 . 23 :

 323 Figure 3.23: Nodes per layer

Figure 3 . 24 :

 324 Figure 3.24: State number ratio between the theoretical worst-case scenario and the Monte Carlo results

3. 7 .

 7 SummaryA Markov Decision Process (MDP)-based path planner for Map-Based Cross-Country Soaring (MBCCS) is suggested in this chapter. The path planner takes updraft maps with uncertainties in the updraft vanish times and updraft strengths as input and calculates planning policies that balance between flight time and the risk of outlandings. A simple flight dynamics model is introduced for propagating the state of the glider. Furthermore, a final glide law is proposed, for calculating airspeed commands during final glide. The technique of Unscented Transform (UT) is utilized in order to transform continuous map uncertainties into uncertain states. The discretization of the outcome probabilities and expected values is done using analytical integration rather than iterative numerical integration in order to reduce the computational load in the decision tree construction. The cost function for the planning problem is based on a glider competition-like rule with penalties for outlandings. The performance of the proposed probabilistic path planner is compared to a state-of-the-art deterministic planner, demonstrating its potential to significantly improve robustness towards updraft map uncertainties thus reducing the risk of outlandings and providing a powerful planning solution for Vision-Based Autonomous Cross-Country Soaring (VBACCS) where updraft maps are inherently uncertain.

2 .FlyFigure 4 . 2 : 2 . 3 .Figure 4 . 3 :

 2422343 Figure 4.2: Flowchart of the blind system architecture (Without remote updraft estimation)

Figure 4 . 4 :

 44 Figure 4.4: Example scenarios with high and low updraft density

  ) where A T represents the size of the weather tile. For the given tile size of 25 km × 7 km, we obtain 65 and 4 updrafts at the defined spacings of 2.5 km and 10 km. During the simulation, vanishing updrafts are replaced by new updrafts, again with random position, lifespan and strength, thus maintaining a constant updraft density over the course of the simulation. Figure 4.4 illustrates the mission setup with varying updraft densities.

Figure 4 .

 4 [START_REF] Bolting | Discrete sliding mode control of small UAS in tight formation flight under information constraints[END_REF] illustrates the cumulated number of outlandings over all Monte Carlo (MC) runs as function of the different image processing noise and updraft density levels. As illustrated in Figure4.5a, at high updraft density, the "blind"

Figure 4 . 5 :

 45 Figure 4.5: Number of outlandings depending on image processing noise and updraft density

Figure 4 .

 4 Figure 4.6 illustrates the three flight trajectories as well as the relevant updrafts. The red solid

Figure 4 . 6 :

 46 Figure 4.6: Flight trajectory comparison (position)

and 4. 7 ,

 7 where both the altitude and airspeed profiles are illustrated. Without remote updraft estimation, the glider's only option in minimizing the distance to the target lies in a direct approach to the target, flying at the airspeed of the best glide ratio, as illustrated in Figure4.7b. Only updrafts that lie on the direct path to the target can be utilized. With no updrafts lying on its

Figure 4 . 7 : 4 :

 474 Figure 4.7: Flight trajectory comparison (altitude and airspeed)

Figure 4 . 8 :

 48 Figure 4.8: Estimated updraft parameters

and 4 .

 4 8h) are additionally illustrated fore the sake of explanation. The gray shaded areas indicate the 3-σ error intervals of the parameters and the dotted vertical lines in all figures indicate the split time. Notice that the updraft estimates are triggered at different times. This is because the vehicle successively perceives the updrafts during the flight. Furthermore, we can see significant variations in the quality of the parameter estimates. These variations result from different observability conditions. Individual clouds are only visible for short periods of time, compared to the entire mission length. This explains why the estimates freeze when no more measurements are available.

  , 4.6b and 4.7a and also in Figures 4.8g

  and 4.8h, by the dashed vertical line. The vehicle has no other option than

and 4 .

 4 8h illustrate that the uncertainties related to the parameters of updraft #3 are of small magnitude which enables a more precise prediction of the action outcome. As a consequence of anticipating the risk in the decision-making, the vehicle accomplishes the mission without an outlanding.

Figure 4 .

 4 Figure 4.9b and Tables 4.4c and 4.4d show the flight time statistics at low updraft density. The blind architecture results in a flight time of 51.73 min for the one successful flight. At both noise

Figure 4 .

 4 Figure 4.10: Average altitude

k 1

 1 weights the cost of a path depending on the remaining distance from the outlanding site to the target. At the beginning of a flight, k 0 will typically have the same impact on all routes, since the destination is far away and not yet reachable. The parameter k 1 induces a more visible impact on the flight trajectories. For illustrating this effect, a flight with equal meteorological parameters was repeated three times, varying the setting of k 1 , i.e. the distance-dependent term, while keeping k 0 at 5000 s. The parameter settings and the corresponding flight trajectories are illustrated in Figure4.11. Note that according to Equation (3.42), the three values for k 1 can be transformed into conveyor belt velocities (Chapter 3) of 0.

Figure 4 . 11 :

 411 Figure 4.11: Trajectories depending on the distance-depending cost function parameter k 1

  three system architectures are established and tested, each of these relying on different updraft estimation and/or path planning methods. Furthermore, four Monte Carlo setups are defined covering two realistic updraft density levels and two image processing noise levels, additionally enabling sensitivity studies of the architectures w.r.t. the most crucial parameters. The Monte Carlo (MC) simulation results reveal the capacity of the proposed remote updraft estimation to significantly reduce the number of outlandings on typical cross-country soaring flights, even when only relying on a deterministic path planner. As expected, the simulation results demonstrate the capability of the proposed probabilistic path planner to further decrease the number of outlandings. Additionally, flight time statistics indicate that the probabilistic path planning approach yields acceptable increases of the average flight time to the destination, compared to the deterministic planner.

F 4 )

 4 where C D , C Y , C L and C l , C m , C n are the dimensionless aerodynamic coefficients for forces and moments respectively. The dynamic pressure is defined as q = 0.5ρV 2 a where ρ specifies the air density. The matrix R b←s conducts transformations between the stability frame and the body frame of the vehicle. Note that the stability frame s is similar to the body-fixed reference frame b. The stability frame is obtained when rotating the body frame about the Y b axis by the angle of attack α. Specifying the rotation rates Ω = [p q r] T and translational velocities V = [u v w] T in the body frame, both angular and translational accelerations are defined asΩ = I -1 (M -Ω × IΩ) = I -1 (M a -Ω × IΩ) , = F a + F g = F a + R b←inAdopting the notation from[START_REF] Beard | Small Unmanned Aircraft Theory and Practice[END_REF], the coefficients for aerodynamic forces and moments as introduced in Section A.1.1 are obtained as the sum of the datum-coefficients, the actuator in-crements, and the body rate damping incrementsC D = C D 0 + (C L 0 + C Lα α) 2 1 πeΛ + C D δe + C Dq c 2V a q, C Y = C Y β β + C Yp b 2V a p + C Yr b 2V a r + C Y δa δ a + C Y δr δ r , C L = C L 0 + C Lα α + C Lq c 2V a q + C L δe δ e , C l = C l 0 + C l β β + C lp b 2V a p + C lr b 2V a r + C l δa δ a + C l δr δ r , C m = C m 0 + C mα α + C mq c 2V a q + C m δe δ e , C n = C n 0 + C n β β + C np b 2V a p + C nr b 2V a r + C n δa δ a + C n δr δ r , (A.5)where e is the Oswald factor and Λ stands for the wing aspect ratio. The actuator deflections for aileron, elevator, and rudder are denoted δ a , δ e , and δ r respectively.The AVL software, as introduced by Mark Drela at the MIT, is used for determining the aerodynamic coefficients of the Cularis glider based on the characteristic geometric properties of the vehicle. Therefore, the geometric properties of the Cularis were measured. The 3D-AVL geometry of the Cularis glider is illustrated in Figure A.1. The identified coefficients are listed in Tables A.2 to A.5.

Figure A. 1 :

 1 Figure A.1: AVL airframe geometry plot of the Cularis glider

L C LA. 1 . 2 .

 12 δa = 0.000 C L δe = 0.773 C L δr = 0.000 Y C Y δa = -0.012 C Y δe = 0.000 C Y δr = 0.187 l C l δa = -0.339 C l δe = 0.000 C l δr = 0.006 m C m δa = 0.000 C m δe = -2.475 C m δr = 0.000 n C n δa = 0.017 C n δe = 0.000 C n δr = -0The Polar Curve As discussed in Chapter 3, the polar curve is a fundamental aerodynamic characteristic for glider planes. It relates the airspeed V a to the sink rate v si of the vehicle and is well approximated by a second order function such that v si = aV 2 a + bV a + c. (A.6) 126 Chapter B-Flight Control Laws p are fed back in order to stabilize the aircraft's roll mode.

2 x2 u, 2 Figure D. 1 :

 221 Figure D.1: Region decomposition in the transformed frame

2 x2λ 1 dx 1 dx 2 c 1 10 ) 12 ) 1 = - λ 1 2 √ 2π k 1 , 1 + λ 2 α 2 u 2λ 1 λ 2 x2 u/l + λ 2 α u β u λ 1 + λ 2 α 2 u,su/l = λ 1 + λ 2 α 2 l 2λ 1 λ 2 λ 1 + λ 2 α 2 l.(D. 15 )

 221101212112121221215 It is now shown how the two integrals c 1 and c 2 can be solved analytically.For c 1 , integration over x1 and rearranging yields Integration over the substituted variables t and s, and rearranging yieldsc u erf tuerf tl -k 1,l (erf (s u ) -erf (s l )) , (D.13)where the corresponding bounds on the new variables aretu/l = λ x2 u/l + λ 2 α l β l λ 1 + λ 2 α 2 l (D.14)and the parameters k 1u/l are defined as k 1,u = e -Note that "erf" specifies the error function which is defined according to erf

u 2 λ 1 +λ 2 α 2 erf αu x2 +βu √ 2λ 1 -erf α l x2 +β l √ 2λ 1 x 2 α 2 u λ 1 +λ 2 α 2 u erf tu -erf tl -α 2 l e - β 2 l 2 λ 1 +λ 2 α 2 l λ 1 +λ 2 α 2 lα u k 5 ,

 2112225 (erf (s u ) -erf (s l )) u erf tuerf tl -α l k 2,l (erf (s u ) -erf (s l )) . (D.[START_REF] Golaz | A Large-Eddy Simulation study of cumulus clouds over land and sensitivity to soil moisture[END_REF]

  gives an overview of the estimator's functionality. It points out the parallel structure of the filter, separating the position estimation from the cloud size parameter estimation. The entry point is the vision block. It provides the pixel positions of the two cloud feature points Y LR for all identified clouds which are transformed into vision measurements Y c . The The vector X d = [d 1 , . . . d l ] specifies a filter bank, i.e. a set of parameter filters that run in parallel. For newly initialized states, the filter bank is extended.Cloud positions X p and cloud parameter states X d are transmitted to the updraft cloud relation block where X u is generated. Note that updrafts for which the parameter estimation has not yet been triggered are not considered in the updraft map. Therefore, the map contains l updrafts. This state reduction is not illustrated in the flowchart for the sake of easier readability.

		position estimation
		prediction	update	augment
			init
	vision	center point	data association	f (X p , X d )
		diameter	init
		prediction	update	augment
	latter are directed to the data association which calculates whether a measurement originates
	from an already initialized (in.), or from an uninitialized cloud (unin.). The correction/update
	step is conducted with measurements from initialized clouds. Uninitialized measurements cloud size parameter estimation
	Y unin.

c,k are used to initialize new cloud positions. The estimated state is augmented by these positions. Therefore, the dimension of X p is time-varying. At the next filter call, state predictions are conducted with the augmented state. The vector Θ contains data association information, necessary to separate between cloud width measurements Y in. d and Y unin.

d .

Table 2 .

 2 .1. They represent a typical low-cost RGB camera.

	Parameter Value	Unit
	x 0	980	[px]
	y 0	640	[px]
	S x	454550 [px/m]
	S y	454550 [px/m]
	f	0.76871	[mm]

1: Camera characteristics

Table 2 .

 2

	2: Process noise parameter values
	Parameter Value Unit
	Σ y LR	2I 4×4 [px]

Table 2 .

 2

	Parameter Value	Unit
	s	5	[-]
	γ	0.085 2 [m/s] 2

3: Measurement noise parameter values

Note that the units refer to the corresponding diagonal elements of the Σ-matrices. In case the same unit applies to several entries, the unit is only noted once. As mentioned earlier, model uncertainties are not considered in this work. Thus, most of the process noise covariance entries are zero. Furthermore, note that the performance of the image processing algorithms, that were suggested in

[START_REF] Halioui | Conception d'algorithme de traitement d'image pour navigation autonome d'un drone planeur[END_REF]

, have not yet been assessed. That said, in this work, the standard deviation of the vision measurement y LR is assumed to be 2 px. An overall system evaluation that also considers increased measurement noise is conducted in Chapter 4.

The remaining filter parameters are noted in Table

2

.4. It specifies the data association threshold for matches s (see Section 2.4.3), and the threshold covariance γ (see Section 2.4.4) used to trigger the parabola tracking.

The sample time of the two filters is selected to be 0.1 s.

Table 2 .

 2

4: Threshold parameter values

Table 2

 2

	.5 lists the initialization

The uncertainty in the initial cloud positions (defined by the values Σ pc xy,i and Σ z cb ) is selected to be 150 m and 10 m respectively. This comparably large uncertainty represents the phase where the vanish time and updraft velocity estimation is not yet triggered. The mean parabola

Table 2 .

 2 5: Initialization parameters for the first Monte Carlo run coefficients e 1 , e 2 and e 3 contained in the vector d are selected to generate clouds with a lifespan of 40 min and a maximum size d max of 750 m on average. Their standard deviation is set to be 10% of the mean initial value. Thus, for each cloud i at run j d ij ∼ N d 0 , Σ 2 d 0 . The vertical velocity scale w is set to 3 m/s and κ is 1/ dmax wu,rz = 0.0033 1/s where wu,rz is 2.455 m/s for a minimum altitude of 500 m and a cloud base altitude of 1200 m.

  .[START_REF] Johnson | A multistation set of whole sky imagers and a preliminary assessment of the emerging data base[END_REF] The first constraint refers to the maximum operational altitude of the glider which we limit by the cloud base altitude z cb . That is, as soon as the glider reaches the cloud base altitude, any climb is interrupted. The second constraint states that the glider does not climb in updrafts with non-positive updraft velocities. Therefore, the exit altitude is inherently equal to or greater than the arrival altitude p z,a and the exit time is equal or greater than the arrival time. The two constraints are summarized asp z,a ≤ p z,e ≤ z cb , t e ≥ t a .For defining a region of possible exit states related to actions a α , the rectangular bounding box (consisting of exit altitude bounds p z,e , p z,ē and exit time bounds t e , t ē) is calculated from the normal distribution s ep ∼ N sep , P se p and some selected confidence interval. The mathe-Figure 3.8: Confidence ellipse of s ep with related bounding box box with the constraints, yielding a saturated bounding box p z,ẽ = min (p z,ē , z cb ) , The region comprised by the saturated bounding box is divided in m × n rectangular bins B, where m and n stand for the number of altitude and time segments respectively. A bin is defined such that B ij . . = {(p z,e , t e ) |p z,e min i ≤ p z,e ≤ p z,emax i , t e min j ≤ t e ≤ t emax The probability of terminating an action from a bin is defined as the normalized integral of the Probability Density Function (PDF) f se p over a region R P s ep ∈ B . . = 1

	p z,e ˜= max p z,e , p z,a , t e ˜= max t e , t a , m×n i=1 R i f se p ds ep	R	f se p ds ep ,	(3.38)
	t ẽ = t ē.
				(3.36)
	saturated bounding box
	confidence ellipse
	sep s a		bounding box confidence ellipse
			se
			s a
	p z,e [m] p z,e [m]		p z,e z cb z cb
			p z,a
	Probabilistic State Transition Model t a		p z,a
			p z,e
	Uncertain estimates of both the vanish time t v and the climb rate w propagate into uncertain t e t e t a t e [s]
	exit states through Equation (3.34). The probabilistic state transition model quantifies out-comes of actions in terms of probabilities and related conditional expectations. t e [s] Figure 3.9: Discretized saturated bounding box

An Unscented Transform (UT) (see Section 2.4.1) is used to propagate the uncertainty P xu into an uncertain exit state s ep and its corresponding covariance P se p through the system equations Equation

(3.34)

. Two constraints are imposed on the exit altitude: matical derivation of the bounding box is shown in Appendix C. It is important to note that since the state s ep features a normal distribution, the ellipse and thus the bounding box can violate the state constraint from Equation (3.36) by comprising regions of the state space that are not accessible. This effect is illustrated in Figure

3

.8. Therefore, we saturate the bounding j },

(3.37)

where 1 ≤ i ≤ m, 1 ≤ j ≤ n and t e min , t emax , p z,e min , p z,emax are the borders of a bin. The discretized saturated bounding box is illustrated in Figure

3

.9 for m = 4 and n = 3 respectively.

  e , t e ) |p z,e ∈ [p z,e min , ∞) , t e ∈ [t e min , t emax ]}, if (p z,emax = z cb , t e min > t a ) {(p z,e , t e ) |p z,e ∈ [p z,e min , ∞) , t e ∈ (-∞, t emax ]}, if (p z,emax = z cb , t e min = t a )

{(p z,e , t e ) |p z,e ∈ [p z,e min , p z,emax ] , t e ∈ (-∞, t emax ]}, if (p z,emax ≤ z cb , p z,e min > p z,a , t e min = t a ) {(p z,e , t e ) |p z,e ∈ (-∞, p z,emax ] , t e ∈ (-∞, t emax ]}, if (p z,e min = p z,a , t e min = t a ) {(p z,e , t e ) |p z,e ∈ (-∞, p z,emax ] , t e ∈ [t e min , t emax ]}, if (p z,e min = p z,a , t e min > t a ) B, otherwise. (3.39)

  In addition to the bin-related probability, we seek to calculate the conditional expectation of a bin. Mathematically, it can be stated asE s ep |s ep ∈ B = 1 P s ep ∈ B R g s ep f se p ds ep (3.40)where the function g saturates the state s ep such that

			bins and constr. regions
			sat. bound. box
			conf. ellipse
			sep
			s a
			
	p z,e [m]	p z,e,sat =	   z cb , if p z,e ≥ z cb    p z,a , if p z,e ≤ p z,a       p z,e , otherwise.	z cb
			
		t e,sat =	  t	p z,a
			a  	t e
			t e [s]
			(a) Unconstrained
			bins and constr. regions
			sat. bound. box
			conf. ellipse
			sep
			s a
	p z,e [m]			z cb
				p z,a
			t a	t e
			t e [s]
		(b) Constrained, follwing Equations (3.38) and (3.39)
		Figure 3.10: Bin-related probabilities

a , if t e ≤ t a t e , otherwise.

Thereby, the saturation function g maps all states s ep that violate the state constraints onto the nearest non-constraint state regions. In Appendix D it is shown how the integral Equation (3.40) can be solved analytically for rectangular regions, thus again being a computationally efficient approach in comparison to an iterative solution calculation.

Table 3 .

 3 .1 (map 1). 1: Updraft map parameter settings A second map with the same mean updraft parameters is generated. For the second map, increased uncertainties are selected by doubling the uncertainties of the first map (see map 2 in

					Updraft
		Parameter Unit #1	#2	#3	#4
		tv	[s]	950 1000 1000 800
		w	[ m s ]	2	2	2	2
	Map 1:	σ tv σ w	[s] [ m s ] 0.15 0.15 0.15 0.15 200 150 100 10
	Map 2:	σ tv σ w	[s] [ m s ]	400 300 0.3 0.3	200 0.3	20 0.3

Table 3 .

 3

[START_REF] Richardson | Upwind dynamic soaring of albatrosses and UAVs[END_REF]

. This enables a parameter study w.r.t. varying map uncertainties. Here, uncertainties apply only to the vanish time and the climb rate while the updraft position uncertainties are supposed to have a negligible effect. The cloud base altitude for both maps is at z cb =1200 m and the minimum allowable altitude is set to z min = 500 m defining an operational altitude range of 700 m.

Table 3 .

 3 .2. 2: Mean and standard deviation of the CTG

	Map Mode Mean [s] Std [s]
	1	det. prob.	8424 5321	19652 46232
	2	det. prob.	33745 8260	15626 17131

  uncertainties, the number of outlandings is reduced from 16 to 8 and from 41 to 16, representing a reduction of 50% and 61% respectively, when relying on the suggested probabilistic path planner.

	Figure 3.18 illustrates the vehicle's trajectories and outlanding sites (noted as a black x) result-
	ing from the two planning policies for low map uncertainties.
		4,000
		2,000		#4		#2	#3
	x [m]	0				#1
		-2,000
		-4,000		Updrafts		Start	Destination	Outlanding site
		0	0.2	0.4	0.6	0.8	1	1.2	1.4	1.6	1.8	2
							y [m]			•10 4
					(a) Deterministic
		-2,000
		0	0.2	0.4	0.6	0.8	1	1.2	1.4	1.6	1.8	2
											•10 4

Table 3 .

 3 3: Monte Carlo settings for generating sample treesIn Figure3.23, both the average number of states per layer (average MC), and the worst-case number of states per layer (worst-case MC), as obtained from a 100-run Monte Carlo study are illustrated in a semilogarithmic plot.

	10 -1	1	2	3	4	5	6

Table 3 .

 3

			TC Runtime [s] BI Runtime [s] Cumulated Runtime [s]
		mean	5.2674	0.3561	5.6235
	MC	std	6.9427	0.0122	6.9548
		wc	25.4021	0.3956	25.7894

4: Runtime statistics

Table 4 .

 4 2 specifies typical meteorological parameters that can be encountered during the summer in the northern hemisphere. The values for both updraft lifespans and strengths, as well as the cloud-base altitude, and the convexing mixing layer thickness z i are taken from[START_REF] Müller | Meteorologische Aspekte des Streckensegelfluges: Thermik von A bis Z[END_REF][START_REF] Allen | Updraft model for development of autonomous soaring uninhabited air vehicles[END_REF].

	The average cloud width amplitude dmax is set to 750 m which according to Equation (2.19)
	converts into a factor κ of 0.0033 1/s. The perturbation gain from Equation (2.6) for generating
	variating updraft strengths is selected to have a standard deviation of 0.12.
	Parameter Distribution	Unit
	z cb	const.,1200	[m]
	z i	const., 2000	[m]
	κ	0.0033	[1/s]
	c	N ∼ (1, 0.12)	[-]
	w	const., 3	[m/s]
	τ u	N ∼ (20, 2)	[min]

Table 4 .

 4 2: Setup-independent meteorological parameter valuesThe major parameter settings for both the estimator and the path planner are listed in Table4.3.Note that the sample time of the remote updraft estimation as well as the path planning parameters are the same as used in Chapters 2 and 3. Also note that the two parameters k 0 and k 1 from Equation (3.42) are selected to be conservative, to prioritize avoiding outlandings over a reduction of flight times.

	Parameter	Description	Value Unit
	t s,est	Remote estimation sampling time	0.1	[s]
	t s,pp	Path planning sampling time	30	[s]
	k 0	Penalty for outlanding (constant)	5000	[s]
	k 1	Penalty for outlanding (remaining distance) 0.01	[-]
	m	Altitude discretization	5	[-]
	n	Time discretization	3	[-]

Table 4 .

 4 3: Major path planning and estimation parameter value settings

Table A

 A Lp = 0.000 C Lq = 10.600 C Lr = 0.000 Y C Yp = 0.016 C Yq = 0.000 C Yr = 0.1411 l C lp = -0.634 C lq = 0.000 C lr = 0.1577 m C mp = 0.000 C mq = -19.446 C mr = 0.000 n C np = -0.058 C nq = 0.000 C nr = -0.042 Table A.3: Rate coefficients Force/Moment w.r.t. δ a w.r.t. δ e w.r.t. δ r

		.2: Datum coefficients
	Force/Moment	w.r.t. p	w.r.t. q	w.r.t. r
	L	C

  e

				-	x2 2 2λ 2	erf	α u x2 + β u √ 2λ 1	-erf	α l x2 + β l √ 2λ 1	dx 2 .	(D.17)
	Using the formula for integration by parts, we can write
	c 2 = -	1 √ π x2,u √ 2 x2,l e √ -2λ 2 λ 2 x2 2	e 2 π √ √	x2 2λ 2 e -(α u x2 + β u ) 2 erf α u x2 + β u √ 2λ 1 2λ 1 α u √ 2λ 1 -erf -e -(α l x2 + β l ) 2 α l x2 + β l √ 2λ 1 2λ 1	x 2,u x 2,l α l √ 2λ 1	dx 2 .	(D.18)

p,k Y in. c,k X u,k X d,k X p,k X d,k|k-1 X p,k|k-1 X d,k X p,kFigure 2.13: Data flow of the estimation process

Acknowledgements

Chapter 3

Path Planning For Cross-Country Soaring

As stated in Chapter 1, path planning is a crucial task for improving the reliability of autonomous cross-country soaring vehicles. However, there has been relatively little research on that field. This is for two reasons. First, autonomous soaring is a young field of research.

Second, to this date, remote updraft estimation functionality does not exist.

The authors of [START_REF] Kahveci | Optimal static soaring of UAVs using vehicle routing with time windows[END_REF] use the Floyd-Warshall algorithm [START_REF] Floyd | Algorithm 97: shortest path[END_REF] to calculate an exact solution of the planning problem. In [START_REF] Kahveci | A heuristic search algorithm for maneuvering of UAVs across dense thermal areas[END_REF], a heuristic based sub-optimal solution of the planning-problem is suggested for long distance flights through areas with high updraft densities. Interestingly enough, in both contributions, full and certain a priori knowledge of all future updrafts over the entire mission horizon are supposed to be available. However, in Vision-Based Autonomous Cross-Country Soaring (VBACCS), the vehicle perceives its environment during flight and the updraft estimation depends on both the vehicle's flight trajectory and the randomness of the nature. Moreover, the remote updraft estimation as introduced in Chapter 2 provides updraft maps of limited size, determined by the estimation performance and the detection range of the image processing device. The mission destination thus may lie far beyond the horizon of the provided map and only a fragment of total flight can be calculated. Most importantly, the updraft maps employed in VBACCS are imperfect, i.e. they feature uncertainties. Ignoring these uncertainties in the planning task can increase the number of mission failures and thus restrict the applicability of VBACCS. These are the main evidences that underline the need of more advanced planning algorithms for VBACCS.

That said, in this chapter, a probabilistic path planner for VBACCS is proposed. Beyond in-Chapter 4

Simulation Evaluation

As pointed out earlier, the main concern of this thesis is to contribute towards more reliable autonomous soaring. Therefore, in Chapters 2 and 3, methods for remote updraft estimation and map-based path planning algorithms are proposed. In the corresponding chapters, both of these contributions are assessed individually. In this chapter, a joint evaluation of the two components is presented. The evaluation is conducted against a variation of crucial parameters. Additionally, the system's performance both in terms of flight time and outlandings is compared to the performance of the state-of-the-art cross-country soaring approach, from the authors of [3].

The estimation and path planning methods presented in this work are platform-independent.

UAV Glider Model

This section introduces the flight dynamics model of a small scale glider plane that is used in Chapters 2 to 4 to generate representative flight trajectories.

A.1. Aircraft Model

Appendix B Flight Control Laws

A set of low level control laws is designed in order to stabilize the Cularis glider from Appendix A. While there exists a variety of advanced controller design methods, here, the purpose of the controllers is simply to generate a realistic glider motion in the simulation and to enable the vehicle to track high-order commands that are generated from both the guidance and the path planning. These commands include the airspeed for the longitudinal motion and the bank angle for the lateral motion.

In the following section, the structure of these controllers is shortly introduced, and the calculated gains are noted.

B.1. Bank Angle and Airspeed Control

The control loops for both bank angle and airspeed were designed using the technique of Eigenstructure Assignment, as explained by the authors of [START_REF] Brockhaus | Flugregelung. Springer[END_REF]. The corresponding control structures are noted in Equations (B.1) and (B.2) respectively.

Note that the controller gains are scaled with the dynamic pressure q. This procedure enables to obtain quasi-equal closed-loop performance across the flight envelope. The reference dynamic pressure qref = 48.33 kg/(ms 2 ) indicates the nominal value at which the linear controllers were designed.

For the bank angle controller, the control error Φ cmd -Φ and its integral, as well as the roll rate

B.3. Thermal Centering Controller

Note that the functionality of the thermal centering controller has been explained in Chapter 4.

Following Equations (4.2) and (4.3), bank angle commands are generated according to

and the commanded turn rate is a function stated as

The steady state radius r ss and the gain K tc specify the control law gains of the thermal centering law.

Appendix C

Bounding Box Calculation

Let x = [x 1 x 2 ] T be a random vector such that

x ∼ N (x, P) ,

with mean x and covariance P. Further, let λ = [λ 1 λ 2 ] T specify the two eigenvalues that are related to the covariance matrix P. Given a scaling factor ξ, as well as the eigenvalues λ, we can specify the semi-minor a el and semi-major b el axis of the related ellipse such that

The scaling factor is related to the selected confidence level 1-α. More specifically, the scale factor ξ is defined by

The rotation angle Θ of the ellipse w.r.t. the coordinate frame is

Then, the bounds of the ellipse's surrounding box are

Appendix D

Bivariate Conditional Expectation

Consider the two dimensional random variable from Appendix C. Further, let f x (x) be the probability density function of the bivariate normal distribution according to Assume the probability P (x ∈ B) is known, for example calculated using the methods provided in [START_REF] Drezner | On the computation of the bivariate normal integral[END_REF]. Then, the double integral in Equation (D.3) has to be solved in order to determine the conditional expectation.

Assume the covariance matrix P can be diagonalized as P = VΛV T , (D.4)