UAS planning and trajectory generation for safe and long-duration oceanic and coastal missions.

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Tesis Doctoral Ingeniería Automática, Electrónica y de Telecomunicación UAS planning and trajectory generation for safe and long-duration oceanic and coastal missions Autor: Leopoldo Rodríguez Salazar Directores: José Antonio Cobano Suárez Aníbal Ollero Baturone Departmento de Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2018 Tesis Doctoral Ingeniería Automática, Electrónica y de Telecomunicación UAS planning and trajectory generation for safe and long-duration oceanic and coastal missions Autor: Leopoldo Rodríguez Salazar Directores: José Antonio Cobano Suárez Profesor (Universidad Pablo de Olavide) Aníbal Ollero Baturone Catedrático Departmento de Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla 2018 Tesis Doctoral: UAS planning and trajectory generation for safe and long-duration oceanic and coastal missions Autor: Leopoldo Rodríguez Salazar Directores: José Antonio Cobano Suárez Aníbal Ollero Baturone El tribunal nombrado para juzgar la Tesis arriba indicada, compuesto por los siguientes doctores: Presidente: Vocales: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Date: Leticia If I lay here would you lie with me and just forget the world Acknowledgements The last three years were very challenging but also rewarding both in my professional and personal life. This thesis is an accomplished goal that could not happen without the support of so many people. Following the example of the Academy Awards, first, I would like to thank my director, Jose Antonio Cobano, who has helped me since day one, not only in aspects related to this thesis but others that have simplified the way through here, he was my mentor and also became my friend. His guidance, expertise, and vision have helped me to tackle different problems and sort all kinds of obstacles in this research and during my stay in Seville. I would also like to thank my thesis director and my tutor Prof. Anibal Ollero Baturone for his guidance, and to push me every day towards excellence and for being a role model of hard work. I want to thank Prof. Tor Arne Johansen for his work in managing the MarineUAS project that has provided us with vast experiences in the field and to organize broad activities that have enriched my professional life. I want to thank my dearest friend Fotios Balampanis who was my partner in crime in this adventure, we were in the same boat, and sailed this rough waters successfully. I hope to share more experiences with you in the future. I would like to give special thanks to Rafael Salmoral, founder of Proskytec; thanks to this project became a terrific friend. He was an integral part of the experiments of this thesis and has helped me enormously with different aspects of the project. I would like to thank Juan Braga, for his support and encouragement on publishing. Margarida Faria, for her great talks and for being a friend of Leticia and me. Thanks to Du Ho, for the great times in all the activities, and thanks to my colleague and friend Bahareh for her friendship and long or not so long conversations. Finally, I would like to thank the rest of MarineUAS fellows, because their companionship, expertise, and ideas helped me enormously in this goal and certainly would help me in the future. You guys are among the most brilliant people that I’ve ever met. My most profound gratitude goes to Leticia, my beloved wife and the love of my life. For her unconditional love and for being my pillar, support and best friend in this adventure, words are not sufficient to thank her for keeping me focused on what matters, for not letting me fall, for telling me when I was wearing the sweater in the wrong way after a long night at work. We have a great future together, thanks for being part of this dream, because my dreams are also yours! Thanks to my mom, who has been my greatest hero and my role model. Thanks again for not letting me miss a single day in elementary school, thanks for teaching me to read and for encouraging me to pursue my dreams. To my grandma that left a few years ago but her example and legacy are still here. Thanks to mom’s and grandma’s cooking lessons Mexico was not that far. And to all the Rodriguez Salazar family, and to the Aguilar Aguilar family, many thanks for your support and for always encouraging me. III XAbstract algorithms; this, together with the chosen software and hardware architectures allows portability to other hardware components and the possibility of incorporating additional functions. Short Contents Agradecimientos V Resumen VII Abstract IX Short Contents XI Notation XVII 1. Introduction 1 1.1. Motivation 1 1.2. Related Work 3 1.3. MarineUAS framework 8 1.4. Thesis contribution and summary 9 2. Wind estimation and identification 13 2.1. Wind models 14 2.2. Ocean influence on the wind field 20 2.3. Equations of motion 26 2.4. Wind vector estimation 29 2.5. Wind feature identification 31 2.6. Wind estimation and feature identification algorithms 37 2.7. Chapter summary and conclusions 43 3. Energy-efficient trajectory generation 45 3.1. Dynamic soaring 45 3.2. Curve selection and parametrization 52 3.3. Chapter summary and conclusions 55 4. Wind efficient flight and trajectory tracking 59 4.1. Area partitioning technique 59 4.2. Guidance and Control Strategies 66 4.3. Moving Frame Analysis 68 4.4. Chapter summary and conclusions 77 5. System integration 79 5.1. Hardware integration 81 5.2. Software integration 84 XI XII Short Contents 5.3. Testing benches 90 5.4. Chapter summary and conclusions 93 6. Systems validation and verification, real experiments and results analysis 95 6.1. Wind field estimation and features identification testing 95 6.2. Trajectory-related simulations 110 6.3. Validation Experiments 120 6.4. Autonomous swooping maneuver experiments 132 6.5. Integration experiments 135 6.6. Chapter summary and conclusions 139 7. Conclusions 141 7.1. Perspectives and Future Work 144 Appendix A. Software versions and hardware characteristics 147 A.1. Software versions 147 A.2. Hardware versions 148 Appendix B. Energy Estimation 151 List of Figures 153 List of Tables 157 List of Codes 159 Bibliography 161 Index 169 Glossary 171 Índice general Agradecimientos V Resumen VII Abstract IX Short Contents XI Notation XVII 1. Introduction 1 1.1. Motivation 1 1.2. Related Work 3 1.2.1. Atmospheric Energy Harvesting 3 1.2.2. Trajectory Generation and Tracking 6 1.3. MarineUAS framework 8 1.4. Thesis contribution and summary 9 2. Wind estimation and identification 13 2.1. Wind models 14 2.1.1. Temporary wind disturbances 17 2.1.1.1. Discrete gust 17 2.1.2. Spatial wind gradient 18 2.2. Ocean influence on the wind field 20 2.2.1. Continuous gust estimation 24 2.2.1.1. Long term prediction 24 2.2.1.2. Short term prediction 25 2.3. Equations of motion 26 2.4. Wind vector estimation 29 2.5. Wind feature identification 31 2.5.1. Weibull Distribution 32 2.5.2. Genetic algorithm 34 2.5.2.1. Fitness function for wind estimation 35 2.5.2.2. Error tests for the fitness function 36 2.6. Wind estimation and feature identification algorithms 37 2.6.1. Wind estimation with direct computation method 37 2.6.2. Wind feature identification algorithm 38 2.7. Chapter summary and conclusions 43 3. Energy-efficient trajectory generation 45 XIII XIV Índice general 3.1. Dynamic soaring 45 3.1.1. The Rayleigh cycle 46 3.1.2. Swooping maneuver characterization 48 3.2. Curve selection and parametrization 52 3.2.1. Curve generation algorithm 54 3.3. Chapter summary and conclusions 55 4. Wind efficient flight and trajectory tracking 59 4.1. Area partitioning technique 59 4.1.1. Waypoint sequencing algorithm 61 4.1.1.1. Intermediate waypoints and control points 64 4.2. Guidance and Control Strategies 66 4.2.1. Trajectory tracking guidance method 67 4.3. Moving Frame Analysis 68 4.4. Chapter summary and conclusions 77 5. System integration 79 5.1. Hardware integration 81 5.1.1. Autopilot and FDES module 82 5.1.2. Hardware communications 83 5.2. Software integration 84 5.2.1. Communication framework implementation 85 5.2.2. Wind field estimation and feature identification 87 5.2.3. Trajectory generation and sequence reconfiguration implementation 88 5.2.3.1. Curve framing implementation 88 5.2.4. Control Implementation 89 5.2.5. Flying modes and control parameters 90 5.3. Testing benches 90 5.3.1. Hardware-In-The-Loop setup 91 5.3.2. Software-In-The-Loop and simulation bench setup 92 5.4. Chapter summary and conclusions 93 6. Systems validation and verification, real experiments and results analysis 95 6.1. Wind field estimation and features identification testing 95 6.1.1. Wind feature identification validation simulations 98 6.1.2. Wind estimation and feature identification SITL tests 108 6.2. Trajectory-related simulations 110 6.2.1. Sequence reconfiguration simulations 111 6.2.2. Trajectory generation and tracking SITL tests 115 6.2.3. Trajectory generation assessment on a full flight plan 116 6.3. Validation Experiments 120 6.3.1. Experimental platforms 120 6.3.2. Experimental sites 122 Índice general XV 6.3.3. Tuning experiments 123 6.4. Autonomous swooping maneuver experiments 132 6.5. Integration experiments 135 6.6. Chapter summary and conclusions 139 7. Conclusions 141 7.1. Perspectives and Future Work 144 Appendix A. Software versions and hardware characteristics 147 A.1. Software versions 147 A.2. Hardware versions 148 Appendix B. Energy Estimation 151 List of Figures 153 List of Tables 157 List of Codes 159 Bibliography 161 Index 169 Glossary 171 Notation αAngle of attack βSideslip angle χcControl law for track angle χkVon Karman’s constant for turbulent zone χfTrack angle ∆Va|zChange of airspeed due to dynamic soaring ηLoad factor εresidual error of fitness function γfFlight path angle γFlight path angle γcControl law for flight path angle κWeibull scaling parameter κcCharnock’s constant for turbulent zone ωpPeak angular frequency ψHeading angle, yaw angle ΦCumulative standard normal distribution φPitch angle ρAir density σ2 nmeasurement noise covariance σu,σv,σwIntensities of the turbulence among the body frame axis θRoll angle τPeriod of a Rayleigh cycle τiTension value for an interpolated spline τL1L1-period υWeibull scaling parameter AAmplitude of the wave AσShear gradient coefficient aTTangential acceleration aCCentrifugal acceleration aAcceleration vector BBinormal vector of Frenet-Serret frame bAccelerometer bias CLLift coefficient CDnDrag coefficients CD0Zero-lift drag coefficients DDrag polar XVII XVIII Notation dmgust length ESystem energy ¯ E(q(t)) Total energy over trajectory e0Oswald efficiency elLateral error in tracking eoLongitudinal error in tracking ewForward error in tracking ggravity constant hAltitude above sea level Hu,Hv,HwDryden transfer functions J(t)wMatrix of temporal wind gradient JwMatrix of spatial wind k(t)Curvature expressed in linear time LLift Lu,Lv,Lwspacial wavelengths LLyapunov function mUAS mass mΩ,nΩCurvilinear coordinates ˜ MiSlope unitary vectors for an interpolated splines NNormal vector of Frenet-Serret frame n Accelerometer measurement noise p,q,rRotational rates in the body frame Q(X,X)Covariance matrix q(X,x)Covariance vector q(x,x)Covariance value SWing area sArc length TThrust TTangent vector of Frenet-Serret frame u,v,wÆirspeed components VaTrue airspeed magnitude VgGround speed magnitude VζSinking speed through the air wWind vector field W20 Reference wind velocity at 20 ft. WλLogarithmic shear effect in the vertical direction WϕPump effect due to the wavelength Wp,Wp(x)Predicted wind speed WuMean velocity profile in the turbulent zone WmGust magnitude Wds minimum wind speed for dynamic soaring Wmp Most probable wind speed at a specific altitude wNED SSteady ambient wind in the NED frame Notation XIX wi(x,y,z)Wind velocity expressed in the inertial frame wshear Horizontal wind gradient ˜ WkFriction velocity in the turbulent zone x,y,zPosition components of the platform z0Transition altitude or reference altitude for wind shear z0,kReference altitude for the friction velocity zTAltitude where the shear velocity is zero 6Chapter 1. Introduction the parameters of the Weibull distribution. Furthermore, there have been techniques for gust detection, which were implemented for UASs: Afridi et al. [ 32 ] have proposed a system for gust detection that aims for the real-time gust detection with a machine learning tool. Yeo et al. [ 33 ] proposed the utilization of onboard flow measurements for gust detection. Etele [ 34 ] has written a complete theory on gust modeling for UAV, which is the basis of a control strategy designed to cope with the unwanted effects of gusts. Wenz et al. [ 26 , 28 , 27 ] and Johansen et al. [ 25 ] have proposed the approximation of continuous gusts with a Dryden wind model, which is used for the estimation of the aerodynamic coefficients. Once the wind is accurately estimated, the conditions necessary for dynamic soaring have to be determined. Bencantel et al. [ 35 ] stated a series of conditions that are necessary for dynamic soaring, which are translated into a minimum shear required, which is proportional to the towing area, air density, and gravity acceleration and is inversely proportional to the aircraft mass. The results show that a more efficient aircraft concerning the aerodynamic parameters requires less shear to perform dynamic soaring. 1.2.2 Trajectory Generation and Tracking Although there have been multiple approaches to solve the guidance and control problem for Unmanned Aerial Systems, there has been few that have proposed strategies for static or dynamic soaring. Nevertheless, few authors have suggested control strategies that deal with this issue. This section discusses the different approaches utilized and shows a mention of the static soaring case since it is a widely studied form of energy harvesting. For thermal soaring, Acevedo et al. [ 36 ] have proposed the utilization of multiple gliders in a distributed approach for thermal detection with a model-based estimation that enables determining what the optimal time in a thermal to extend the duration of the flight is. Cobano et al. [ 10 ] have proposed the utilization of a trajectory planner based on the RRT* (Optimal Rapidly-exploring Random Trees) for accurate thermal exploitation. Clarke and Chen [ 37 ] evaluated different methods for trajectory generation for thermal exploitation, and have recommended the use of the Piggott [ 38 ] and the Reichmann [ 39 ] techniques for thermal centering. The Reichmann’s method can be summarized into three premises: 1) as climb improves, then the flying circle shall be flattened. 2) as climbing deteriorates, the circle must be steeped, and 3) if the climb is constant, keep constant bank. These techniques have been widely used for manned gliders and were also taken into consideration for autonomous vehicles [ 40 ], in which a control system was designed to achieve the behavior mentioned above. 1.2 Related Work 7 In the case of dynamic soaring, several research efforts deal with the planning and the trajectory generation. Lawrance [ 41 ] presented a method that utilizes Gaussian Process Regression (GPR) to generate a wind map that works together with a path planning algorithm to generate a trajectory that aims for energy gain. The trajectory considers an optimal solution without actually looking for a specific type of soaring but for maximizing a reward function that contains information of the wind field. Gao [ 11 ] has stated techniques for both thermal exploitation and shear exploitation in the context of surveillance missions for UAS’s. Gao has proposed the execution of a dynamic soaring cycle, i.e., traveling across the shear layer of the wind. To merge the cycles, Gao presented a Dubins-path- based trajectory planner. Akhtar [ 42 ] addressed the problem of generating problems for energy exploitation with dynamic soaring with a direct method in real time, in which the trajectories are dependent on the coefficients of a polynomial that define the paths. The process proves to have advantages (mainly those related to real-time implementation) over the classic approach presented by Zhao [ 43 ], who consider the dynamic soaring problem a nonlinear optimal control problem and has generated off-line solutions for different wind gradient profiles. Zhu [ 44 ] suggested the determination of a so-called equilibrium position, in which the potential energy gained while flying throughout the vertical wind gradient, is higher than the wasted kinetic energy of flying upwards. Syama et al. [ 45 ] stated an off-line theory, in which given a suitable wind profile, specific characteristics of the aircraft are incorporated into a dynamic optimization routine that utilizes numerical methods to determine the optimal trajectory. Bower [ 46 ] defined a trajectory optimization method for dynamic soaring, with the objective of finding a solution for optimal periodical trajectories in different wind fields, which are stored in the database. However, Bower focused part of his work on determining the platform characteristics which are optimal for dynamic soaring. This work was also considering large-area oceanic missions, which required the study of the interaction between the sea and the ocean. The trajectory planning and control could be executed in a multi-stage approach, for instance, with the utilization of a receding horizon control to track such trajectory. In low level, there is an inner-controller that deflects the control surfaces. The details on the implementation of the receding horizon control can be found in [ 47 ]. Montella and Spetezer [ 15 ] proposed a reinforcement-learning approach that controls the aircraft while flying in shear conditions. Again this was a multi-stage approach with a predefined trajectory that is flown with a sample-based controller. The selected method for trajectory tracking was proposed by Miao et al. [ 48 ], in which the generated curve is framed with a moving frame, e.g., the Frenet-Serret frame or the Type 1 bishop frame. After the curve framing process, the error is determined and is used to define a Lyapunov function that monitors the positional and angular errors between the 8Chapter 1. Introduction body frame of the vehicle and the moving frame of the curve. Finally, this function is utilized to provide with flight path angle commands and also with track angle commands. 1.3 MarineUAS framework The work presented in this thesis is part of the MarineUAS project, which is a program part of the Marie Sklodowska-Curie actions of the European Commission to strategically strengthen research training on Autonomous Unmanned Aerial Systems for Marine and Coastal Monitoring. It is a comprehensive researcher training program across a range of industrial and education partners in several countries, including, Norway, Spain, Portugal, Sweden and the Czech Republic. The program is designed to have a high impact on the training of individual Early Stage Researchers (ESR’s) and their knowledge, skills and their future careers. MarineUAS has established a unique cooperative environment. It takes benefit of the partners’ extensive and complementary expertise, operational field experience, and experimental facilities [49]. Figure 1.1 Research and training activities on the MarineUAS project include the realization of tests in harsh environments, training on RPAS piloting, aiming to design the most innovative solutions for the problems of coastal and oceanic monitoring. The research lines of the project include the design of fault-tolerant control for coastal monitoring, iceberg detection, performance optimization, model-based diagnosis, cooperative motion planning, among others. The project aims to build a solid foundation for long-term European excellence and innovation in this field by sharing research infrastructures for field testing and disseminating the research and training outcomes and best practice of MarineUAS into the doctoral schools of the partners, as well as by fostering long-term partnerships and collaboration. 1.4 Thesis contribution and summary 9 The motivation of this research is aligned with the motivation of the project, that considers the complexity of monitoring the coasts and Exclusive Economic Zones (EEZ) of the European Union, that stated in the European Strategy for Marine and Maritime Research, that there is a need to protect and manage the vulnerable natural environment and marine resources in a sustainable manner. The conditions for monitoring these resources are sometimes harsh and demanding. Therefore, innovative ideas and approaches, including the use of UASs, are necessary to deal with this objective. 1.4 Thesis contribution and summary This thesis presents the development and implementation of a system capable of extending the flight duration of UASs by utilizing the available energy from the atmosphere by two means: the first one is an off-board wind-based waypoint sequencing planning or an on-board win-based waypoint sequencing reconfiguration. The second is the selection of portions of the given flight plan, that is suitable to execute dynamic soaring trajectories in a safe manner. Compared to the reviewed methods in Section 1.2.2, the system focuses on the mission and optimizes portions of it to gain energy and reduce the fuel/battery consumption by harvesting energy from the atmosphere. Also, there was a large experiments campaign, which included flights in a variety of locations, conditions and scenarios. The objective was to validate the presented methods at different stages. One critical aspect is that the system was developed with a systems-engineering approach, in which a set of requirements is defined and cascaded to every function, which permits to determine clear means for validation and verification of the different system functions. The goal was to prioritize the use of COTS components and the real-time implementation of the selected methods, which drive most of the design decisions of the research. These are the main contributions of this thesis: The development of a fully functional system for the implementation of smart AEH in a surveillance mission context. The combination and improvement of different methods for wind estimation and feature identification that are normally utilized in other fields. A novel methodology for trajectory generation with parametric curves that are suitable for dynamic soaring. The integration of a novel adaptive control algorithm for tracking complex curves that allow the realization of dynamic soaring. 10 Chapter 1. Introduction The consideration of a simpler reconfiguration of the waypoints sequencing as a means to achieve more efficiency combined with the dynamic soaring. This thesis consists of seven chapters, with complementary information in two appendices. Next, an overview of the contents of each chapter is presented. Chapter 2 explains the wind models utilized for the wind field estimation and the characterization of features of interest (shear and gusts), together with the sea-air interaction; it also explains in detail the methodology for determining the wind relevant characteristics that are necessary to dynamic soaring. Chapter 3 describes the methodology employed for trajectory generation, the bioinspiration for the selected type of trajectories, the parametrization technique and the trajectory generation algorithm. Chapter 4 sets out the methodology for path planning and trajectory tracking. It starts from the methodology to create or reconfigure a wind efficient flight path, followed by the methodology for the selection of the control points for the parametrized trajectory and finalizing with the proposed tracking methodology, including a definition of a tracking error definition and the control methods. Chapter 5 explains the integration of the different methods described in the previous chapters into a functional architecture. Moreover, the software and hardware architectures are described, together with the architecture of testing benches that were utilized for validation and verification activities. Chapter 6 shows a large experimental campaign, that was performed throughout the development of the preceding system, including the results and the analysis from the Hardware- and Software-In-The-Loop experiments, as well as the field experiments. Chapter 7 has the conclusions about the research and presents directions for further research. Several publications in international conferences and indexed journals have been generated in the course of this thesis. These publications are: Leopoldo Rodriguez, Jose A. Cobano, Anibal Ollero Wind characterization and mapping using fixed-wing small unmanned aerial systems, 2016 International Conference on Unmanned Aircraft Systems (ICUAS), pp. 178-184, Arlington (VA, USA), June 7-10, 2016. 1.4 Thesis contribution and summary 11 Leopoldo Rodriguez, Jose A. Cobano, Anibal Ollero, Wind field estimation and identification having shear wind and discrete gusts features with a small UAS, 2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 5638- 5644, Deajeon (Korea), October 9-14, 2016. Leopoldo Rodriguez, Jose A. Cobano, Anibal Ollero, Small UAS-Based Wind Feature Identification System, Part 1: Integration and Validation Sensor, Special issue on UAV- Based Remote Sensing, December, 2016. Leopoldo Rodriguez, Jose A. Cobano, Anibal Ollero, Qualitative Analysis of Guidance and Control Methods for Wind Gradients Exploitation with Small Fixed Wing UAS. The 2017 International Workshop on Research, Education and Development on Unmanned Aerial Systems, RED-UAS 2017, 102-107, Linköping (Sweden), October 3-5, 2017. Leopoldo Rodriguez, Jose A. Cobano, Anibal Ollero, Architecture of a Flight Endurance Enhancement System for Maritime Operations with Fixed Wing UAS, ROBOT 2017: Third Iberian Robotics Conference, pp. 171-182, Seville (Spain), November 22-24, 2017. Leopoldo Rodriguez, Jose A. Cobano, Anibal Ollero, UAS Soaring Trajectories Considering the Atmospheric Boundary Layer for Oceanic Long Duration Missions, 2017 MTS/IEEE OES OCEANS conference, Anchorage (Alaska, USA), 2017. Leopoldo Rodriguez, Fotios Balampanis, Jose A. Cobano, Ivan Maza, Anibal Ollero, Energy-Efficient Trajectory Generation with Spline Curves Considering Environmental and Dynamic Constraints for Small UAS, 2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 1739-1745, Vancouver (BC, Canada), September 24–28, 2017. Leopoldo Rodriguez, Jose A. Cobano, Anibal Ollero, Smooth Trajectory Generation for Wind Field Exploitation With a Small UAS, 2017 International Conference on Unmanned Aircraft Systems (ICUAS), pp. 1241-1249, Miami (FL, USA), June 13-16, 2017. Leopoldo Rodriguez, Fotios Balampanis, Jose A. Cobano, Ivan Maza, Anibal Ollero, Wind Efficient Path Planning and Reconfiguration of UAS in Future ATM, Twelfth US- A/Europe Air Traffic Management Research and Development Seminar (ATM2017), Seattle (Washington, USA), 2017. 12 Chapter 1. Introduction Fotios Balampanis, Leopoldo Rodrıguez, Jose A. Cobano, Ivan Maza, Anibal Ollero, A path enhancement experimental framework for heterogeneous non holonomic flying vehicles Submitted to IROS 2018 - Under Review. Leopoldo Rodriguez, Jose A. Cobano, Alberto Galán, Guillerno Heredia, Anibal Ollero Wind estimation and trajectory generation methods for UAS long duration missions, Submitted to Journal of Intelligent Robotic Systems (JINT 2018) - under review 2 Wind estimation and identification The wind of change blows straight into the face of time. Scorpions Kn owledge of the environment is critical for UAS mission planning and execution. In particular, the atmosphere plays a vital role while flying; information about the wind phenomena could help to improve the safety of the mission and even to a possible increase of the flight duration. Hence, different methods to understand the wind interaction with the UAV shall be studied to provide with information to avoid any hazardous condition and to increase (to take advantage of) the available energy. This chapter describes a multi-stage approach for the characterization of the wind dynamics, starting from the identification of the incident wind field on the aircraft while flying followed by a characterization of the features that allow the creation of a wind map for features of interest that can be exploited for energy gain. As stated in Chapter 1, there have been numerous research efforts on estimating the wind vector of a moving aerial vehicle. The most common approach is the utilization of low-cost navigation sensors, such as the combination of IMU, GNSS, and an air probe through the use of a Bayesian filter or even directly with a simplified model. For this research, a simple approach namely, direct computation method is selected as the means for estimating the wind vector. Furthermore, this method is compared with existing and popular approaches which have been proposed not only for soaring applications but others, such as icing detection. Also, the knowledge of the wind field, i.e., wind velocity, the rate of change of wind velocity and the wind gradient, is of particular interest to perform safer and more efficient flights. The additional information on the wind may allow better navigation and to identify potential areas to harvest or to avoid. This identification process starts with the local wind vector estimator; then this information can be stored and utilized to create a full wind map with identified features of interest that can be used online or off-line for different purposes. 13 14 Chapter 2. Wind estimation and identification Surface Layer Planetary Boundary Layer Free atmosphere region Convective mixed layer Entrainment zone Ground Figure 2.1 Structure of the Planetary Boundary Layer (PBL) [11]. 2.1 Wind models The atmosphere is divided into five layers based on temperature: troposphere, stratosphere, mesosphere, thermosphere, and exosphere. The phenomena studied in this research occurs in the troposphere, which contains up to 75% of the air on the planet. This layer extends from the surface up to an altitude of that goes between 9 km and 16 km of height, depending if it is measured at the poles or the equator respectively. The troposphere can be divided in the Planetary Boundary Layer (PBL), also called Atmospheric Boundary Layer (ABL), and the free atmospheric region. According to Gao [ 11 ], the thickness of the PBL varies with temperature, so that in extreme cold it can go up to 100 m and in very hot equatorial areas can be of several kilometers. Also, the PBL can be divided into three layers where different phenomena occur: the surface layer, the convective mixed layer, and the entrainment zone (see Figure 2.1). There is a significant presence of wind gradients in the surface layers, which are determined by different factors. These include surface friction, or interaction with the ocean waves, daytime solar heating, radiative night cooling, humidity, stability conditions, etc. This zone and the interaction with the surface and the ocean is where most of the research on this chapter is focused because these conditions affect the flight characteristics of small UAS. Conjointly, the mixed layer is also of interest due to the presence of convective wind patterns, namely thermals or thermal columns, which are suitable for energy harvesting. The forthcoming models are expressed in different reference frames : 1. The Inertial frame, I : it has its origin at the center of the earth. The x axis points to the vernal equinox, the z axis points to the north polar axis, and y points to the resultant right-handed direction. 2.1 Wind models 15 Prime meridian To Vernal Equinox Equator ye yI xI xe zIze N E D xB yB zB Local level tangent plane Figure 2.2 Different reference frames utilized in the definition of the wind models. 2. The Earth-Centered-Earth-Fixed (ECEF) frame, e : it has its origin at the center of the earth. The x axis points to the intersection of the prime meridian with the equator, the z axis points to the north polar axis and y points to the resultant right-handed direction. 3. The navigation frame, or North-East-Down (NED) frame: it has its origin at the aircraft center of mass. Its axes are oriented along the geodetic directions defined by the Earth’s surface: one axis points to the polar direction. Another axis points east along a latitude curve, and the last axis points towards the earth surface. 4. The body frame B : it has its origin at the platform UAV center of mass, the x axis points to the right wing, the zaxis points towards the bottom of the vehicle 5. The wind frame W : it is a convenient frame to express the aerodynamic forces and moments acting on an aircraft, having its origin on the aircraft aerodynamic center with the x axis pointing forward and aligned to the velocity vector. The positive y axis passes through the right wing and the axis points downward. The aircraft flow angles (angle of attack and sideslip) define the relationship between the wind and the body axis system. These frames are illustrated in Figure 2.2 but the wind frame. Note that the wind frame is illustrated in section 2.3 Langelaan, Alley, and Neidhofer, in [ 18 ], proposed an approach that estimates the wind vector together with the gradients and acceleration. This calculation was considering the 22 Chapter 2. Wind estimation and identification ments performed at a altitude of z0=19.5 m . The spectrum equation as a function of the angular frequency ω: SPM(ω) = 8.1×10−3g2 ω5exp −5 4ωp ω4!(2.14) The peak angular frequency ωp can be expressed in terms of the measurements at an altitude of 19.5 meters: ωp=2πfp≈0.8772 g W19.5 (2.15) The amplitude, A, of the wave in the sea surface is derived. A=r2Z∞ 0 SPM(ω)dw ≈7.4×10−2W19.5 g(2.16) It is assumed that the ocean behavior is that of a fully developed sea, i. e. steady wind over long period with waves considered as a single sinusoidal wave. Then, the wave surface height Z(x,t)is Z(x,t) = Acoskx +ωpt(2.17) in which: k=ω2 p g(2.18) Since the wind velocity is dependent on the shape of the wave, an additional reference frame is needed. In such surface-fitted Ω , the curvilinear coordinates (mΩ,nΩ) are defined as mΩ=x+Aekz sinkx −ωpt(2.19a) nΩ=−z−Aekz coskx −ωpt(2.19b) 2.2 Ocean influence on the wind field 23 The partial derivatives of (mΩ,nΩ) are necessary to add them up to the horizontal and wave gradient logarithmic influences ∂mΩ ∂x=1+Akekz coskx −ωpt(2.20a) ∂nΩ ∂z=Aκekz sinkx −ωpt(2.20b) (2.9) should transformed to be expressed with Ω frame coordinates so that the shear profile follows the contour of the wave and considering a logarithmic profile, i.e. 1<A<2 , Wshear =˜ Wk χlogn z0(2.21) The 3D wind velocity and the wind spatial gradient can be derived from the previous relationships. The wind velocity is expressed in terms of the logarithmic shear effect in the vertical direction Wλ , and in terms of the pump effect due to the wave Wϕ . First, considering (2.21) the wind velocities with the added up shear effect are: Wxλ=W(u)∂nΩ ∂x(2.22a) Wzλ=W(u)∂nΩ ∂z(2.22b) The velocity perturbations due to the wave pumping are defined as: Wxϕ=−A√gκcosκx−ωpte−κn(2.23a) Wzϕ=−A√gκcosκx−ωpte−κn(2.23b) Resulting in the following 3D wind equations for the boundary layer expressed in the wind frame: Wx=Wxλ+Wxϕ(2.24a) Wy=0(2.24b) Wz=Wxλ+Wxϕ(2.24c) 24 Chapter 2. Wind estimation and identification The spatial gradient models shown in (2.8) and the particular case for the oceanic gradient model shown in (2.21) may allow the extrapolation of the measurements to other altitudes. Typically, this is performed based on an empirical power law considering an exponential distribution between two wind measurements at a certain altitude. However, the exponent is subject to the region stability and the roughness, which is addressed in the shear model. In the presence of temporary gradients, i.e. discrete or continuous gusts, the distribution gets wider, however in the long term, it still may be able to approach to that measurement. 2.2.1 Continuous gust estimation Two methods are proposed for continuous gust identification in both short and long term. Due to the nature of the continuous gust phenomena, there is high variability in the conditions, and not always are suitable for soaring. Therefore, the objective of the continuous gust estimation is to provide with information to determine if the gust is suitable for soaring or to avoid the area. The considered methods should be suitable for real-time implementation. Also, the computational cost have to be relatively low for an on-board implementation. 2.2.1.1 Long term prediction The first method utilizes a standard Gaussian Process (GP) regression [ 11 ] to provide a suitable estimation in the long term. The purpose of this estimation is to determine if there is a potential gust. Considering a set of vertical wind observations of size ¯ M , ˆ Wz=ˆ Wz,i|ˆ M i=1 . The wind speed prediction ¯ Wp(x)at any location xcan be expressed as: ¯ Wp(x) = ¯ M ∑ i=1 kiˆ Wz,i(2.25) in which ki is the i -th coefficient of the linear combination of wind measurements ˆ Wz . Based on the work presented by Park et. al. [ 57 ], an optimal coefficient is determined by minimizing the prediction error. min kE(¯ Wp(x)−Wp(x))2=min k(kTQ(X,X)+σ2 nIk−2kTq(X,x)+q(x,x)) (2.26) which can be determined by calculating the covariance matrix Q(X,X) and the covariance vector q(X,x) between every two observations at locations X and x ; finally q(x,x) 2.2 Ocean influence on the wind field 25 represents the covariance value . This leads to express standard GP regression of the linear predictor as: ¯p(X) = kTˆ Wz=q(x,X)Q(X,X)+σ2 nI−1ˆ Wz(2.27) and the covariance value Q(¯p(x))) can be expressed as: Q(¯p(x)) = q(x,x)−q(x,X)Q(X,X)+ σ2 nI−1q(x,X)(2.28) where σ2 nis the measurement noise covariance. 2.2.1.2 Short term prediction The second approach is used to perform long-term predictions with a non homogeneous regression prediction model. Lerch and Thoraninsdottir [ 58 ] have performed a comparison between three non-homogeneous regression model, which allows to produce predictions in day time window. Since the intention of the intended testing flight campaigns is to store data into a single database, a big amount of data can be utilized to perform the predictions with the selected regression model. At this stage, the truncated normal model was selected as a form of wind estimation. Being W the wind speed and X1,...,Xj the ensemble member forecasts, the predicted distribution of Wcan be approximated by a truncated normal distribution: W|X1,...,Xj∼N[0,∞)(µ,σ2)(2.29) where the mean µ is an affine function of the ensemble forecast and the variance σ2 is an affine function of the ensemble variance. If these exchange members are exchangeable [ 58 ], the distribution function of the Truncated Normal (TN) distribution F(z) is given by: F(z) = Φµ σ−1Φz−µ σ(2.30) for z>0, where Φis the cumulative standard normal distribution. This is indeed a simple non-homogeneous method. However, results indicate that the training period to produce accurate predictions in one-day ahead forecasts is of the equivalent 30 days of continuous measurements [ 58 ]. Hence, the goal is to be able to forecasts gusts at least after a month of experimental campaign, which intend to have a very bast amount of data for the prediction models. In addition the long term prediction model goes for measurements that are accumulated in the period of 6 months [ 58 ], from which analysis is considered a goal for the end of the experimental campaign. These two 26 Chapter 2. Wind estimation and identification methods will allow to understand the gust formation in a region and also to characterize the region stability which will allow to have a better wind forecast for the soaring applications. After the understanding of the estimation models, the next step is to have a set of equations of motion that will be utilized in the simulation stage and also to understand certain aspects that are critical to the trajectory generation methods. 2.3 Equations of motion The understanding of both spatial and temporary gradients are the base for the estimation of the wind field and the identification of features of interest for possible soaring. First, a set of equations of motion that consider the spatial gradient need to be defined, for this, the research work presented by Bowler [ 46 ] compared two approaches of equations of motion, one by Sachs [ 59 ] and the other by Zhao [ 43 ], both in the context of dynamic soaring. The Sachs model has fewer trigonometric terms and it does not require the wind rate to be solved, however it presents a disadvantage when it comes to set the initial conditions in a trajectory optimization problems. In contrast, Zhao’s approach considers the use of airspeed, flight path angle and heading angle as state variables, which allows a more intuitive solution. As in [ 46 ], the set of equations of motion that will be utilized is the one presented by Zhao. Consider the aircraft depicted in Fig. 2.6 located at point r with coordinates (x,y,z) in the inertial frame; the acceleration of a point mass can be derived considering acting lift L , drag D , thrust T , and weight mg forces. The wind field velocity components (Wx,Wy,Wz) are defined in the inertial frame. The angle γa corresponds to the air,-relative climb and ψais the heading angle of this air-relative motion vector. First, the vehicle velocities are defined in the inertial frame: ˙ r=va+w(2.31) the components of ˙ rare defined as: ˙x=Vacosγcosψ+Wx(x,y,z,t)(2.32a) ˙y=Vacosγcosψ+Wy(x,y,z,t)(2.32b) ˙z=−Vasinγ+Wz(x,y,z,t)(2.32c) 2.3 Equations of motion 27 R xi yi zi ^ ^ ^ mg L D xi ^ yi ^ Va -zi ^ ψa γa Wz -Wy -Wx T Figure 2.6 Relationship between the inertial and the wind frame, illustrating the most significant angles and the forces acting on the platform. The differentiation with respect of time of the velocities in (2.32a) - (2.32c) are the body axis accelerations accelerations: ¨x=˙ Vacosγcosψ−˙ γVasinγcosψ−˙ ψVacosγsinψ+˙ Wx(2.33a) ¨y=˙ Vacosγsinψ−˙ γVasinγcosψ−˙ ψVacosγcosψ+˙ Wy(2.33b) ¨z=−˙ Vasinγ−˙ γcosγ+˙ Wz(2.33c) The lift L and drag D in Figure 2.6 are defined as functions of the true airspeed and the aerodynamic coefficients: L=1 2ρV2 aSCL(2.34) where ρis the air density, Sis the reference area and CLis the lift coefficient. 28 Chapter 2. Wind estimation and identification The drag polar can be defined as an n -degree polynomial, however, the approach proposed in [ 46 ] and [ 43 ] is selected due to their previous use in dynamic soaring calculations: D=1 2ρV2 aS(CD0+CD1CL+CD2C2 L+CD3C3 L+CD4C4 L)(2.35) in which CDn,n=1...3 , correspond to the drag coefficients, which can be approximated with a wind tunnel characterization and also could be estimated with a CFD analysis. With the lift and drag definitions, one can apply the second Newton law to define the sum of forces in the three axis of the inertial frame: m¨x= (T−D)cosγcosψ−L(sinγcosψcosφ+sinψsinφ)(2.36a) m¨y= (T−D)cosγsinψ−L(sinγcosψcosφ−cosψsinφ)(2.36b) m¨z=−(T−D)sinγ−L(cosγsinφ)+mg (2.36c) where φis the roll angle. Now, the wind velocity rate of change in equations (2.33a) - (2.33c) correspond to the spatial and temporal wind gradient. The spatial gradient can be defined in its matrix form Jwas: Jw=    ∂Wxλ ∂x+∂Wxϕ ∂x0∂Wxλ ∂z+∂Wxϕ ∂z 000 ∂Wzλ ∂x+Wzϕ ∂x0∂Wzλ ∂z+Wzϕ ∂z     (2.37) The temporal gradient, which correspond to the stochastic portion defined in (2.3) , is defined as J(t)w=    ∂Wxλ ∂t+∂Wxϕ ∂t 0 ∂Wzλ ∂t+∂Wzϕ ∂t     (2.38) The equations of motion, with the aircraft position at point R , can be then derived in terms of the gradients 2.4 Wind vector estimation 29 L=m cosφ    Va dγa dt +gcosγa−   cosψasinγa sinψasinγa cosγa    T (Jw˙ R+J(t)w)    (2.39a) ˙ γ=1 Va     L mcosφ−gcosγ+   cosψsinγ sinψsinγ cosγ    T (Jw˙ R+J(t)w)    (2.39b) D=1 2ρV2 aSCD(α)+ L2 1 2ρV2 aSπRe (2.39c) ˙ Va=−D m−gsinγ−   cosψcosγ sinψcosγ −sinγ    T (Jw˙ R+J(t)w)(2.39d) ˙ ψ=1 Vacosγ      sinψ −cosψ 0    T (Jw˙ R+J(t)w)    (2.39e) 2.4 Wind vector estimation An initial estimation of the local wind can be provided as in [ 18 ] from the set of equations of motion (2.39) , in which the wind velocity is computed directly from the sensor information, this is referred as the direct computation method. This comes from the relationship of the wind velocity observed in (2.32a) - (2.32c) . If the wind terms are solved, it can be obtained: Wx=˙y−Vacosγcosψ(2.40a) Wy=˙y−Vacosγcosψ(2.40b) Wz=˙z+Vasinγ(2.40c) The airspeed components can be renamed as: u=Vacosγcosψ , v=Vacosγcosψ and w=Vasinγ . Considering the information provided by an IMU, which corresponds to the body axis accelerations (¨xb,¨yb,¨zb), the angular rates (p,q,r), and the attitude angles pitch θ , roll φ , and yaw ψ , and considering a and accelerometer bias b and with measurement 30 Chapter 2. Wind estimation and identification noise n the accelerations considering the wind components expressed in the body frame and the airspeed components in the body frame would be: ¨xb=˙ Wb x,I−˙ub I+qw −rv +gsinθ+bx+nx(2.41a) ¨yb=˙ Wb y,I−˙vb I+ru −pw −gcosθsinφ+by+ny(2.41b) ¨zb=˙ Wb z,I+˙wb I+pv −qu −gcosθcosφ+bz+nz(2.41c) Since the calculation of the rate of change of the velocity, with respect the air mass, and expressed in the body frame, cannot be determined with the on-board sensors it can be estimated with a second order numerical differentiation. Therefore, the wind speed rate of change at the previous step k−1results in    ˙wb x ˙wb y ˙wb z   k−1 =   ¨xb−bx ¨yb−by ¨zb−bz   k−1 +   −gsinθ gcosθsinφ gcosθcosφ   k−1 −   qw −rv pw −ru qu −pv   k−1 −1 2∆t   uk−uk−2 vk−vk−2 wk−wk−2    (2.42) The estimation of the wind vector using the information of the IMU and the groundspeed components from the GNSS, allow a direct computation of the vector. The main disadvantage of this estimation method is that the airspeed components expressed in the body frame are necessary to perform the computation. Most of the COTS autopilots provide an approximation of these values. By looking at the main source code of open license low cost autopilots, one can infer that the estimation is poor, as it only considers two axis and not the full rotations. This approximation was introduced by Premerlani [ 30 ], in which a body-to-earth DCM is defined with a simple derivation, the wind speed measurements based on the known speed triangle. The source code of the selected autopilot has such implementation as seen in Code 2.1. The open source project allows the possibility of implementing the aforementioned method in the selected platform, which is discussed in chapter 5. Code 2.1 PX4 Wind Calculation Code. if (_tas_data_ready && (_imu_sample_delayed.time_us - _airspeed_sample_delayed.time_us < (uint64_t)5e5)) { _state.wind_vel(0) = _state.vel(0) - _airspeed_sample_delayed. true_airspeed * cosf(euler_yaw); 2.5 Wind feature identification 31 _state.wind_vel(1) = _state.vel(1) - _airspeed_sample_delayed. true_airspeed * sinf(euler_yaw); } else { // If we don’t have an airspeed measurement, then assume the wind is zero _state.wind_vel(0) = 0.0f; _state.wind_vel(1) = 0.0f; } An alternative method is proposed by Wenz [ 26 , 28 , 27 ] as part of the MarineUAS project, and Johansen [ 25 ], in which an optimal estimator is used for the estimation of wind velocity as a part of an algorithm that aims to estimate the angle of attack α , the sideslip angle β aerodynamic coefficients of the aircraft while flying. The most popular approach is the utilization of an EKF, which drives to some challenges, such as the problem of distinguishing the changes of the aerodynamic coefficients, for instance in an icing episode, and the changes on the wind velocity. The method presented by Wenz evolved from an EKF that needed persistence of attitude excitation for convergence, to a MHE, the architecture of the system adopted in this research, that contains MHE is shown in Figure 2.7, illustrating the interaction for different functionalities. The initial hypothesis is that the values obtained with the direct computation method, computed with the integration of (2.42) , are sufficiently accurate to provide with the necessary information to identify the features that are suitable for energy extraction. 2.5 Wind feature identification The behavior of the wind is a complex phenomena, which is influenced by different factors that may produce sudden changes that imposes restriction on the modeling. These factors include the terrain (surface stability and surface roughness), the time of the day, time of the year, among others. However, the potential identification of features, such as gusts and the characterization of the shear, will allow the determination of portions of the mission that can be utilized for dynamic soaring. This section describes the algorithms that were proposed and utilized for the identification and characterization of these features, which allow the creation of a live wind map. One of the objectives of the determination of these algorithms is to understand the influence of the ocean in the determination of the adequate conditions for AEH. 38 Chapter 2. Wind estimation and identification values are received, they are processed to obtain a wind estimation that will be published to be stored in the temporary and permanent database. Algorithm 1 Wind field estimation with the direct computation method 1: procedure WindEstimationDCM(ˆa,ˆr,ˆu,ˆ b,↑n,Attitude, RotationRates, t) 2: while INS and GNSS and Pitot is OK do 3: [ax,ay,az]= GetAccelerations(); 4: [˙x,˙y,˙z]= GetGroundVelocities(); 5: [u,v,w]= GetAirspeedComponents(); 6: [θ,ρ,ψ]= GetAttitudeAngles(); .These angles can be obtained as in (2.40). 7: [p,q,r]= GetRotationRates(); 8: [u,v,w]k−2= GetPrevRotationRates(); .Determine at k−2as in (2.42) 9: [wx,wy,wz]= CalculateWind(); 10: [˙wx,˙wy,˙wz]k−1= CalculateWindRates(); .The rates are calculated at step k−1. 11: PublishWind(); 12: PublishWindRates(); 13: end while 14: end procedure Algorithm 1 depends on the measurements computed by the autopilot, or any external source that calculate the functions (line 11 and 12 in Algorithm 1). The communication of the different variables has to allow the real time implementation. The actual communication means is discussed in Chapter 5. The solution of the navigation is calculated in the autopilot with the use of a Kalman Filter. Therefore, the error at each measurement is contained within the autopilot, and since the relationship of the variables is direct, the same error is propagated to the wind estimate. This means that the error will not grow and if the autopilot provides reliable estimates, then the error in the wind velocity vector estimation is contained. Nevertheless, there are checks of the covariance, which can be extracted on the autopilot as part of the looping conditions of the measurement. 2.6.2 Wind feature identification algorithm The wind feature identification algorithm is a prediction tool that aims to detect and identify features that are necessary for the energy extraction process. The algorithm performs a statistical test to a group of wind velocity estimates to determine if a feature is present. First, the module requests a wind estimation and stores it into a temporary database. If there are sufficient estimates (at least 30) from the current flight, the system starts a feature detection process by ordering the wind database with respect the UAS altitude. Since the altitude reading varies with time, even in small amounts, the estimates 2.6 Wind estimation and feature identification algorithms 39 are grouped according to a reference altitude and they are selected if they are sufficiently close to the reference value. Normally the references altitudes are integer numbers and the groups are conformed by those readings between a tolerance of ±1 m. Then, the module should find the corresponding Weibull parameters for each altitude using GA from the grouped wind estimates. If the system finds the Weibull parameters, a most probable wind speed at the corresponding reference altitude is generated. The process is repeated until the local maximum altitude is reached. At first, the system performs an analysis to determine the presence of a shear, which is a very common feature [ 31 ]. The wind prediction module tries to characterize the coefficient Aσ in (2.9) , that minimizes the error between the most probable wind speeds for the reference altitudes. If the estimates are distributed according to the Weibull distribution and there is a value of Aσ that produces an acceptable error into the system (during the testing, the Aσ coefficient was selected when the average error among the different altitudes was ε≤5 m/s), an alert is triggered and the system recognizes the presence of a shear. Afterwards, the system performs a statistical analysis to determine anomalies (significant jumps) in consecutive wind estimates. These were performed by looking for sudden increases into the running standard deviation of the wind estimates. If there is a sudden increase an initial alert is generated because potentially a gust is identified. If the system is not capable to determine accurately the Weibull parameters of the system, most probable wind speeds cannot be fitted into a shear, and/or the running standard deviation presents drastic changes, i.e. there are continuous increases in the running standard deviation, the system assumes the presence of a continuous gust which triggers a short term Gaussian Regression process to characterize the feature. Algorithm 2 describes the insights of the wind feature detection and identification algorithm. Algorithm 3 defines the procedure of gust detection, Algorithm 4 describes the process of grouping altitudes for statistical tests, and finally Algorithm 5 shows the steps of the Weibull parameter estimation with the use of a GA. The goal is to avoid a big computational impact with the implemented GA algorithm (see Figure 2.9). In this case, the GA implementation is part of the function that calculates the most probable wind speed at a certain altitude described in Algorithm 5, the parameters were carefully selected, based on Liu’s approach [ 31 ], which prove a low computational impact, since only the shaping parameter κ is estimated. The other parameters necessary to calculate the most probable wind speed are relationships that are not dependent on the GA, which is a very good choice for the real time implementation of the algorithm. 40 Chapter 2. Wind estimation and identification Algorithm 2 Wind feature detection and identification algorithm 1: procedure RequestWindEstimation 2: WindEst=RequestWindEstimate() .Request from Algorithm 1. 3: WDb = StoreinDB(WindEst) 4: end procedure 5: procedure DetectFeature(WDb) .If database has at least 30 elements 6: if WDb.Size≥30 then 7: Start = True; .Start detection of features. 8: else 9: Alert=InsufficienElements 10: end if 11: if Start=True then 12: WDb = OrderAltitudes(WDb); .Order WDb based on altitude. 13: for i=1←AltMax do .Check for altitudes 1 m to maximum altitude . 14: NearAlts=FindNearAltitudes(WDb,i,thres); 15: AdNearAlts=RequestDb(i); .Additional WDb elements. 16: NearAlts=[NearAlts:AdNearAlts]; .Group elements. 17: WindVelMP=FindMPWVel(NearAlts) .most probable wind speed at i. 18: MPS(i)=Store(WindVelMP); .Store the most probable wind speeds. 19: end for 20: Aσ=CalcCoeff(MPS) .Calculate Aσcoefficient from (2.9). 21: if A.Exist = True then 22: Alert=ShearDetected; 23: end if 24: i++; 25: if Exist(A)=False then 26: GustDetected = GustDetection() 27: end if 28: if GustDetected = True then 29: end if 30: end if 31: end procedure There are two functionalities in Algorithm 2 that need to be defined, one is the gust 2.6 Wind estimation and feature identification algorithms 41 detection algorithm GustDetection(), and the algorithm that groups the altitudes based on an altitude reference FindNearAltitudes(), that are shown in Algorithm 3 and Algorithm 4 respectively. Algorithm 3 Gust Detection algorithm 1: procedure GustDetection 2: DetectJumps(WDb.Velocity,Std(WDb.Velocity)) 3: if JumpDetected=True then 4: JumpCounter++; 5: end if 6: if JumpCounter≥threshold then 7: Alert=ContGustDetected 8: else 9: Alert=DiscGustDetected; .Gust Detected. 10: end if 11: if CommHandler.Request.Alert.DiscGustDetected=True then 12: Gust=DetectJumps.Jumpsize 13: else if CommHandler.Request.Alert.DiscGustDetected then 14: ContGust=PerformGaussianRegressionWDb 15: end if 16: end procedure The algorithm that is used to group the altitudes based on a reference is shown in algorithm 4: Algorithm 4 Grouping near altitudes algorithm 1: procedure FindNearAltitudes(WDb,alt,thres) . Find altitudes in WDb close to alt. 2: Counter = 0; 3: for i=1←WDb.Size do 4: if alt-thres≤WDb(i).Altitude≤alt+thres then 5: NearAlts(Counter++)=WDb(i); .Store whole WDb. 6: end if 7: end forreturn NearAlts; 8: end procedure 42 Chapter 2. Wind estimation and identification The determination of the most probable wind speed at a given altitude is described in algorithm 5 in which ,. Algorithm 5 Weibull parameter calculation algorithm 1: procedure FindMPWVel(NearAlts) .Find altitudes. 2: κ=CalcKappa(NearAlts) .Calculate shaping parameter using GA. 3: ν=1 Mean(NearAlts) Γ(1+1 κ).Calculate scaling parameter from (2.46a). 4: end procedure 5: procedure CalcKappa(Altitudes) .GA Implementation. 6: PopulationSize=50; 7: FunctionTolerance=1×10−3; 8: MaxGenerations=100; 9: CrossOverFraction=0.8; 10: StdAlt =Std(NearAlts); .Calculate standard deviation. 11: MeanAlt=Mean(NearAlts); .Calculate mean. 12: PopKappa==rand(PopulationSize); .Initialize with random population. 13: while ε>FunctionTolerance do 14: for j=1←PopKappa.Size do 15: Results(j)=ObjFunc(PopKappa(j),StdAlt,MeanAlt); .Evaluate Objetive Function. 16: end for 17: Parents=Selection(Results,PopKappa); .Selection of elements for newGeneration. 18: Reproduction(Parents,PopKappa,MaxGenerations;) .Creation of new population. 19: Crossover(CrossOverFraction); .Scattered crossover function. 20: Migration(); .Gaussian Mutation function. 21: end while 22: end procedure The aforementioned algorithms are sufficient to provide the necessary functionalities for the detection and identification of features of interest for the dynamic soaring. The detection process can be executed in real time with the presented algorithms, however, they are based on the simplest implementation and do not consider entirely the effects of the interaction of the atmospheric/ocean boundary layer. Nevertheless, this effects may 2.7 Chapter summary and conclusions 43 not impact significantly in the detection process, but they are relevant for the positioning of auxiliary waypoints that will define the soaring trajectories in oceanic operations. The determination of the wave amplitude and frequency on fully developed seas, based on the wind measurements, will allow the determination of the places where the platform may sink optimally for a higher energy gain. The process of trajectory generation and tracking is discussed in the following chapters. 2.7 Chapter summary and conclusions This chapter presented the different models that were utilized to find a method that extends the flight duration of a UAS. These models were subject of a thorough analysis and were evolved from simple models widely used in aeronautic applications (CS-25/FAR- 25 Airworthiness Standard and some Military specifications), to more complex model that consider the interaction of the ocean with the atmosphere. This evolution is key to understand the importance of having accurate wind measurements and how the different phenomena can be utilized correctly to get the maximum energy gain. Also, the impact of the surface accidents, such as mountains, buildings, etc, is discussed, because this factor will determine the ability to extend the wind predictions to wider areas. In addition, this chapter presented the dynamic model that will be utilized along this work. This model gives big importance to the wind estimation, allowing a full consideration of the impact of the wind estimations into the model. The most important part of this chapter is the estimation methods, from which the direct computation method was selected to be implemented due to its simplicity and the proven results in literature. However, other models are discussed, which were developed as part of other research activities of the MarineUAS model that can be implemented to improve the estimation framework. Furthermore, the methods and algorithms for identification and characterization of features algorithms are presented, including the implementation of the different statistical tests, based on the Weibull distribution. The identification of the Weibull parameters given a dataset was accomplished with the use of a GA and the insights of its implementation was also discussed. The chosen algorithms provide a solid base for the accurate estimation, but also to a real-time implementation which is key to the process of extracting energy. The main contribution of this chapter is the utilization of the Weibull distribution in real-time applications, because it was utilized commonly for very long datasets for wind turbine placement. In addition, the models that consider the interaction of the wind and the ocean are not widely used, specially in the prediction process, which is also an advantage in the selected approach because not only wind predictions can be performed but also the 44 Chapter 2. Wind estimation and identification wave behavior could be inferred from the wind measurements which would be a great advance for autonomous soaring applications. The following chapter will discuss the next stage on the development of a framework for flight duration extension: the trajectory generation in which the different steps for generating a wind efficient trajectory are discussed. 3 Energy-efficient trajectory generation Don’t quack like a duck, soar like an eagle. Ken Blanchard Th is chapter describes the trajectory generation methods that were employed for the energy-efficient trajectory. Firstly, an overview of the different forms of dynamic soaring is presented, with an analysis of the albatross flight, including an energy analysis followed by the details on the implementation of a generator of dynamic soaring trajectories in real time focused in oceanic operations using curve parametrization. The goal of the trajectory generation functionality is to create on a waypoint-to-waypoint basis trajectories that will aim to maximize the efficiency. 3.1 Dynamic soaring The dynamic soaring technique consists in the exploitation of spatial gradients in a wind field, allowing savings in the energy that is utilized to generate lift. This technique can be observed in nature with the albatross (Diomedea Exulans). The phenomena was firstly studied in the late XIX and early XX centuries, with the main purpose to replicate its extraordinary ability to soar very long distances in the Southern Ocean without wind flapping [ 19 , 59 , 17 , 23 ]. The understanding of the dynamic soaring and the possible large-scale replication has its basis on the observation of the albatross (see Figure 3.1). This birds are biologically designed to fly long distances with a remarkable wing span that goes up to 3.7 m on males and a weight that goes from 7 kg to 10 kg [55]. According to [ 23 ], two theories aim to explain the long-distance soaring of the albatross. Although both theories suggest that the albatross exploit spatial gradients, the source of this gradient differs from one to the other. The first one suggests that the albatross utilizes the wind shear, whereas the other indicates that the albatross uses the updrafts formed with 45 46 Chapter 3. Energy-efficient trajectory generation Figure 3.1 The wandering albatross in flight. For a complete illustration of the dynamic soaring an explanatory video can be found in [64] . the interaction in the air/ocean boundary layer. The wind in the air/ocean boundary layer has a wind shear component and a vertical motion, so the contributions of each effect are hard to be studied. Also, the gusts, lee eddies and rolls (temporary gradients) affect the study of each effect since these phenomena have often wind shear and updrafts. Typically, the study of the dynamics of the albatross has to be simplified [ 59 , 22 ]; the most common assumption is that the winds are purely horizontal and the absence of waves, therefore, the wind-slope soaring component is excluded. Moreover, the typical wind models consider at most a logarithmic or an exponential gradient model, which does not recognize the wave effects on the gradient and makes it difficult to understand the flight dynamics of the albatross fully. Richardson [ 23 ] proposes a model based on the observation of the albatross to estimate the change of energy on the albatross flight and also to create a more accurate model of the wind-shear soaring. This model is the combination of the observations made by Rayleigh [ 19 ] that introduced the concept of a bird soaring on the spatial gradient. Also, it is based on the observations of Pennycuick [ 21 ], that observe in detail the soaring process of the albatross in the lee-side of the waves, and also in the dynamic model introduced by Lissaman [ 22 ] that creates one of the firsts dynamic models applicable to both birds and vehicles but with the limitations mentioned above. 3.1.1 The Rayleigh cycle The flight pattern that describes the trajectories of the albatross when performing energy extraction is often referred as the Rayleigh cycle [ 23 ]. This cycle consists in the utilization of two wind layers and it is often modeled based on the premise that a minimum wind speed is required to perform this maneuvers. The aerodynamic forces studied in Chapter 2, are dependent on the airspeed, however there are studies [ 22 ] that offer a comparative of the calculations of the energy while using airspeed or ground speed. There is an increase of airspeed when dynamic soaring is performed, and this increase occurs with the wind- 3.1 Dynamic soaring 47 shear layer cross. However, there is a ground speed increase, which is uncorrelated to the airspeed, that occurs with a heading turn from headwind to downwind. With these statements there is a slight complication in the interpretation of the energy conservation principle. The most popular dynamic equations that offer a solution to this problem were stated by Sachs [ 59 ]. The proposed model has been used widely [ 55 , 46 ] while trying to study and imitate the soaring trajectories. The results show that the albatross needs to get sufficient energy to compensate the drag losses, and this energy comes from both wind-shear soaring and wave-slope soaring. From the observations performed int he southern hemisphere, concretely in Cape Town Southafrica, Richardson [ 23 ] studied albatross soaring in a wind field with gusts that go from 5 m/s to 10 m/s . The soaring pattern, that lasts in average 10 s was described as follows: Each swoop began with a fast flight parallel to and just above the windward side of a wave. This was followed by a turn into the wind and climb of around 10 m - 15 m followed by a downwind descent towards another wave and a turn parallel to the wave. [23] Richardson developed an schematics of this pattern based on his observations (see Figure 3.2), which is a simplification of the real pattern. For simplification, the waves are shown uniformly, illustrating a fully developed sea, in which the input of energy to the waves from the local wind is in balance with the transfer of energy among the different wave components. Considering the shear model presented in (2.9) , it can be observed that with a logarithmic growth ( 1<A<2 ), the majority of the increase in wind speed occurs in the boundary layer. There are measurements that suggest that the critical wind speed increase is located within 2 m over the sea surface. As it can be seen in Figure 2.1, the structure of this area has additional lee eddy effects and gusts. These lee eddies are utilized by the albatross to perform the dynamic soaring, by entering and leaving the boundary layer. The Lee eddy zone is illustrated in Figure 3.3 Several authors [ 21 , 59 , 55 ] refer the lee eddy dynamic soaring as gust soaring. However, this is not purely gust soaring since the energy gain occurs by transitioning between the lee-eddy zone to over the boundary layer. In the actual gust soaring the gusts are aligned with the lift vector togenerate an energy gain. Another form of dynamic soaring, often referred as hill or slope soaring occurs due to the interaction of the wind and the terrain. The updrafts are produced as the wind circulates over hills and mountain ridges [ 45 ]; the actual soaring maneuvers are skewed helical close to the leeward side of the ground. The 54 Chapter 3. Energy-efficient trajectory generation As stated before, four Bézier control points are defined in terms of the curve parameters. The initial and final control point correspond to the consecutive respective waypoints of the path Riand Ri+1: P0=Ri(3.15) P1=Ri+Sin(ξi,ηi+1) 3τi||Ri+1−Ri||ˆmi(3.16) P2=Ri+1(3.17) P3=Ri+1 Sout(ξi,ηi+1) 3τi+1||Ri+1−Ri||ˆmi+1(3.18) The proposed polynomial for a cubic interpolation of the control points is defined as: R(t) = K3t3+K2t2+K1t+K0;(3.19) where the Knparameters are defined as functions of the control points Pn: K0=P0(3.20a) K1=3∗(P1−P0)(3.20b) K2=3∗(P0−2P1+P2)(3.20c) K3=−P0+3(P1−P2)+P3(3.20d) An example of a generated curve of four waypoints can be observed in Figure 3.7. 3.2.1 Curve generation algorithm The calculated curves have to be generated on a waypoint-to-waypoint basis allowing the reshaping of the trajectory in accordance with the wind conditions to maximize the energy gain (see Figure 3.7). The generation algorithm employed in the overall trajectory generation and tracking solution considers both the bank angle method, i.e. determining the bank angle for dynamic soaring and perform semicircular trajectories within a given period and diameter for a desired soaring speed, and the usage of hobby splines for complexer cases. Furthermore, the placement of the set of mission waypoints for energy maximization will be studied in the following chapters. In the case of semi-circular Rayleigh cycle Section 3.1.2, the process is straight forward: if the winds allow the possibility of performing stable 3.3 Chapter summary and conclusions 55 2.5 Distance (m x 103) 2 1.5 1 0.5 0 -1 1 -0.5 0 Distance (m x 103) 0.5 1 0 0.5 1 Distance (m x 102) 2 3 4 Figure 3.7 Example of parametrized curve showing the effects of modifying parameters. The curve is a closed one, showing lower tension and velocity values in black and stiffer parameters in red. dynamic soaring, 10knot ≤W≤30knot , then the period is calculated together with the optimal radius for the soaring speed that does not exceed the maximum load factor. If the winds are higher than the ones mentioned above, then between two consecutive waypoints, one or two additional waypoints will be introduced based on the wind speeds and the distance. This additional waypoints will be the transition between the climbing and descent phases of the Rayleigh cycle and will be positioned based on the desired period which can be estimated initially by determining a desired increment on airspeed based on (3.4) . However, in the case of a fully developed ocean, the period of the Rayleigh cycle will match the wave period determined in (2.15) . The curve velocity functions as in (3.11) are determined so that the maximum curvature does not exceed the semi-circular trajectory which will have the maximum airspeed gain; the tensions τi at each waypoint is chosen to be close to 1 for closed Rayleigh cycles and close to 0 for open Rayleigh cycles. The unitary vectors are chosen to follow the desired flight direction between Rayleigh Cycles. A Rayleigh cycle with a parametrized curve can be observed in Figure 3.8. 3.3 Chapter summary and conclusions This chapter presented a biomimetic approach to replicate the flight of the albatross toextract energy. Firstly an analysis of the albatross flight is presented. These observations provide with different types of trajectory that will result into a different behavior of 56 Chapter 3. Energy-efficient trajectory generation wind direction X Position Y Position Z Position waves flight direction Figure 3.8 Example of a Rayleigh cycle following the method of parameterizing the curves. The desired trajectory is marked in black. The additional waypoint is observed in the middle of the curve and the overall vehicle trajectory follows the dotted flight direction, which is accomplished by selecting the unitary vectors parallel to the desired direction. In this case the tension values were selected to be τ=0.04 for every waypoint . the energy gain process. In addition, the necessary conditions for dynamic soaring are discussed toprovide with sufficient information for the trajectory generation. These trajectories, which are referenced as Raleigh Cycles, are complex curves which are difficult to replicate. Hence, a method that generates a parametric curve that may replicate this trajectory is presented. In other research activities, the dynamic soaring trajectories were considered circles or portions of a circle, which is not always efficient, specially in oceanic conditions. Therefore, the parametric curve offers versatility to produce complex, but smooth curves which really imitate the albatross flight. The parametrization method, offer numerous advantages, because the determination of the control points depends on the desired characteristics of the curve, which could be constrained by the kinematic characteristics of the platform. Moreover, the need to add auxiliary waypoints into the construction of a Rayleigh cycle was identified and address with an algorithm that permits their determination considering the wind information extracted in the previous chapters, that considers the interaction between the ocean and the atmosphere. Therefore, the waypoints are placed in such way that a portion or a full Rayleigh cycle is constructed but it prioritizes the preexisting waypoints, which is key in the context of coastal and maritime surveillance missions. The trajectory generation functionality is the key part of the energy extraction process 3.3 Chapter summary and conclusions 57 of the dynamic soaring, and the generation method is a key contribution of this work. The smoothness, which is a factor that depends on the curvature, is such that it is ensured that the tracking of these trajectories will not represent a problem in terms of flyability. The next chapter discusses the trajectory tracking and the waypoint sequencing reconfiguration, which are the two steps that are performed tomaximize the flight duration without compromising the mission integrity, which most of the time depend on the precise fly-over of the mission waypoints. 4 Wind efficient flight and trajectory tracking If you have everything under control, you’re not moving fast enough. Mario Andretti Th is chapter presents the algorithm employed to perform changes in an initial waypoint sequence to maximize the wind energy gain in a given mission. In the context of a surveillance mission, the typical approach is to generate a set of waypoints, considering a partition technique to cover a specific area. The waypoints can be positioned based on different approaches to optimize the path in terms of distance, probability of finding a feature of interest, etc. An additional step could be performed to set up the sequence of waypoint in the most efficient manner considering the wind. The implementation of a re-sequencing algorithm will permit energy gains in straight flight, and to implement the soaring cycles in better conditions by allowing simpler Rayleigh cycles for higher energy gains. Also, it will prevent the execution of Rayleigh cycles in portions of the flight path that will add unnecessary complexity to the mission by determining if the energy gain will be significant. Initially, a brief analysis on the selected method for complex area coverage is presented, from which the flight plan is generated. Following by the waypoint re-sequencing algorithm. Finally, a methodology for trajectory tracking is proposed and analyzed to complement both the trajectory generation presented in Chapter 3, and the waypoint re-sequencing function. 4.1 Area partitioning technique The trajectory tracking technique was based on the approach presented by Balampanis et al. in [ 67 ], in which a framework for coverage considering complex areas is presented. The 59 60 Chapter 4. Wind efficient flight and trajectory tracking Figure 4.1 The trapezoid represents the field of view of a mounted sensor in the platform, the triangle is utilized as the partition unit of the area, which ensures that at slight changes of attitude as in (a), which change the trapezoid shape, the information will not be lost. This work as presented by Balampanis et. al. [ 67 ] as part of the MarineUAS project.. complex areas could include holes (restrictions) and could be in any shape. Besides, the selected methodology allows the partitioning considering multiple vehicles, which could also be utilized to extend the wind feature identification methods presented in chapter 2. In this research, the Constrained Delauney Triangulation (CDT) [ 68 ] is utilized as the area partitioning technique, and the Lloyd optimization is used to further improvement of the resulting triangles by improving the angle between cells. The cells are chosen to be triangular because with such form it is assured that each cell can be contained into a sensor footprint, as in any surveillance mission, the typical pattern of a sensor footprint is trapezoidal as it is observed in Figure 4.1. The CDT introduces edge constraints to create a mesh of triangles. The circumcenter of each triangle is considered a waypoint to provide with full coverage of the area. A cost is given to each triangle to sequence the flight plan, and this cost can be based on different criteria (task, area or agent-based). In the research presented by Balampanis et al. in [ 67 ], the waypoint planning problem is a graph search problem of some nodes organized in the CDT. The advantage of this approach is that the region can include complex area restrictions, resulting in even more complex areas and the coverage of the area is still sufficient. For the purpose of this research project, which includes the surveillance of complex maritime/coastal areas, the studied regions considered air traffic restrictions in different locations. For instance, regions such as the great Seattle area (see Figure 4.2), were studied due to the amount of restrictions (urban areas, islands, military operations areas, etc.). These regions were partitioned using CDT. The size of the sensor footprint is variable, and depends on the UAS attitude and altitude; this has a noticeable impact on the cell size. For standardization purposes and also following the current regulations on 4.1 Area partitioning technique 61 Figure 4.2 A test case area, north of Seattle. The red polygon defines the mission area whereas several other restrictions apply. Remote-Piloted Aircraft System (RPAS)’s, the cells are typically calculated at an altitude between 100 m and 300 m. 4.1.1 Waypoint sequencing algorithm After the triangulation, the areas are subdivided based on the number of UAS’s. The initial weighting criteria for each cell is the proximity to the border which may result into a spiral path. Furthermore, additional weighting can be provided, such as the one presented by Balampanis et al. in [ 69 ], in which the weighting technique is improved for a variation of a pure pursuit tracking method. This was achieved by weighting the resultant angle between any two waypoints which, will then maximize or minimize smoother curves. However, to adjust this method to the requirements of this research, a reconfiguration technique is proposed in which, once the wind features are identified, the sequence is reconfigured with two main objectives: To maximize the number of legs in which the platform flies with tail-wind, meaning that there is a significant positive component between the predominant wind direction and the leg direction. To identify the potential legs in which a Rayleigh cycle can be performed for energy maximization. Regarding the first objective, there might be a slight complication with the initial weighting problem: if the main weighting contributor is the comparison between the track and the wind direction, then there might be legs with sharper turns. The second objective 62 Chapter 4. Wind efficient flight and trajectory tracking Wind > 10 Knots Figure 4.3 An example of leg reconfigurations where the blue lines represent the sequence generated with the CDT algorithm with a normal weighting, and the red lines sequence represent the sequence reconfigured. The solid red lines represent normal leg trajectories, and the dotted lines represent trajectories that will be flown with a Rayleigh cycles since there is head wind. involves an intricate process; in the Rayleigh cycle description in subsection 3.1.1, there are different cycles that can be performed according to the scenario. There are cycles that may result into a dramatic speed increase, which may produce additional difficulties in the control strategy or in the structural integrity of the vehicle. On the other hand, there are cycles that may reduce the speed even with a significant energy gain. Therefore, an algorithm that considers these restrictions allowing firstly a smart waypoint sequencing, and a selective trajectory generation is needed. The generated curves of the Rayleigh cycles may introduce additional complexity because the energy gain depends very much on how accurate is the tracking therefore, any strategy for creating complex trajectory has to prevent any overload in the control strategy. A graphical example of the re-sequencing algorithm can be observed in Figure 4.3. The wind-based waypoint sequencing reconfiguration occurs after two events are completed: the first is that there is a valid waypoint sequence, and the second is that there is a sufficient number of wind estimates or sufficient wind information. The availability of wind information imposes the first restriction, that is that the sequencing needs to be performed after the vehicle is airborne, since the main information comes from the wind velocity estimates, however, there are still ways to determine wind direction during planning. The determination of when the waypoint list can be reconfigured into a wind-based sequence t0 , occurs when the predominant wind speed and direction is identified at the 4.1 Area partitioning technique 63 mission altitude. This can occur by two means, either with the extrapolation of meteorological data from a near meteorological station. This information can be extracted from the National Oceanic and Atmospheric Administration (NOAA), the European Climate Assessment and Dataset (ECAD), which is very useful in case of oceanic missions, and also from the Meteorological Terminal Aviation Routine Weather Report (METAR) reports available. These lasts are valid if the flight is close to a station that emits these reports. With the initial wind direction, the weighting of the cells in the CDT can be adjusted, giving more weight to the neighboring cells that produce a leg that will have a positive component. The weight can be increased or decreased proportionally to the angular difference between the leg and the predominant wind direction. The predominant wind direction can be adjusted on-board, allowing dynamic changes of the sequence. The process of sequence reconfiguration is illustrated in Algorithm 6. Algorithm 6 Wind Feature Estimation 1: procedure ReconfigureSequence(W,R).Depends on the wind and waypoint sequence. 2: if Wis valid then .Wind information is valid. 3: if Ris valid then .Valid waypoint sequence 4: NewSeq = R 5: for i=0to Size(R)do 6: N=GetNeighbor(R(i)); .Get closest waypoints to the current 7: for j=0to Size(N)do 8: A.ang=CalculateAngles; .Calculate leg angles of neighbors. 9: for k=0to Size(A)do A.weight(k) = A.ang(k)/W.dir .Weight the angles according to the wind direction 10: if A.weight(k)>Thres then 11: M=CheckAngles(A(k)); .Check angles of neighbor waypoints 12: M=Max(M) .get the maximum angle of the neighbor waypoints 13: NewSeq(i+1=R(Index(M))); .adjust waypoint sequence 14: WptReconfFlag(i)=True; 15: end if 16: end for 17: end for 18: end for 19: end if 20: end if 21: end procedure The implementation of the reconfiguration process makes a comparison of different 70 Chapter 4. Wind efficient flight and trajectory tracking κ=dT ds 1 ˙ R(t) (4.15) τ=−1 ˙ R(t) dB dt ·N(4.16) The curve presented in Figure 4.5 is framed utilizing the Frenet-Serret moving frame, the result is shown in Figure 4.6. Y Position (m) 2000 1000 -1000 0 4 1800 1600 1400 3 1200 X Position (m) 1000 800 2 600 400 1 200 0 -200 100 0 200 300 400 600 500 Z position (m) Curve Tangent Normal Binormal Figure 4.6 The curve presented in Figure 4.5 showing the tangential T , normal N and Binormal B vectors after framing it with the Frenet-Serret frame among the different waypoint segments. To move to the Bishop frame, the angle θF should be determined by using the following equation: θF(t) = Zτ(t)˙ R(t)dt +θ0(4.17) Then, the two curvature parameters of the bishop frame can be determined from (4.7) : k1=κ q1+tan2θf (4.18) k2=tanθFk1(4.19) The generated curve, framed on the type 1 bishop frame, is shown in Figure 4.7. The next step is to define and calculate the error between the body frame B and the curve frame F . The error definition, the construction of the Lyapunov function, and the definition of the control laws are based on the work presented in [48]. 4.3 Moving Frame Analysis 71 1800 1600 1400 1200 1000 4 800 X position (m) 3 600 400 200 2 0 -200 1 -2000 0 Y Position (m) 100 0 600 500 400 300 200 2000 Z Position (m) Curve Tangent N1 N2 Figure 4.7 The curve presented in Figure 4.5 showing the tangential T , and the normal N1 and N2 vectors after framing it with the type 1 Bishop frame among the different waypoint segments. T N1 N2 xb yb zb xi yi zi R O P r p e Figure 4.8 Represented error among the curve, the blue frame represents the inertial frame i with origin point at O . The black frame represents the body b frame of the platform with the center of gravity located at R . The red frame represents the frame located at the tracking point P . The e distance represents the positional error between the tracking point and the current aircraft position at R. Consider the frames depicted in Figure 4.8. The goal is to determine the error equations that describe the positional error, e , between the vehicle position at time, t , and a tracking point, P , that will be changing with the aircraft move and also the angular errors between the body and the frame. In the case of Figure 4.8, the curve frame F is a type 1 Bishop frame, showing the N1,N2and tangential Tvectors. 72 Chapter 4. Wind efficient flight and trajectory tracking Consider a rotation matrix CF I that changes a coordinate set from the I to F frame. Such matrix can be parametrized locally by the angles χf and γf that represent the attitude at the tracking point (track and flight path angle respectively). This matrix is: CF I(χf,γf) =    cosγfcos χfcosγfsin χf−sinγf −sin χfcos χf0 sinγfcos χfsinγfsin χfcosγf   (4.20) The rotational rate at point P , ωf , expressed with respect the inertial frame is defined as ωF=   0 ˙ γf 0   +   cosγf0−sinγf 0 1 0 sinγf0 cosγf      0 0 ˙ χf   (4.21) Given a parametrized curve P(t) , in the 3-dimensional space, with the equation as in (3.19) , the arc length, s, can be defined as a function of the coordinates of a tracking point P,(xf,yf,zf): ds =sdxf dt 2 +dyf dt 2 +dzf dt 2 dt (4.22) Then, the values of ˙ γf and ˙ χf can be expressed as a function of the arc length ˙s of the parametric curve, and ˙ χf=χf ds ds dt = darctan ˙xf ˙yf ds ˙s=˙xf¨yf−¨xf˙yf ˙x2 f+˙y2 f3 2 (4.23a) ˙ γf=γf ds ds dt = darctan ˙zf q˙y2 f+˙x2 f ds ˙s=¨zf˙x2 f+˙y2 f−˙zf¨yf˙yf+¨xf˙xf q˙y2 f+˙x2 f˙x2 f+˙y2 f+˙z2˙s(4.23b) The velocity at the tracking point Pcan also be expressed in the curve frame Fas dp dt F =   ˙s 0 0   (4.24) The velocity at the point R, expressed in the inertial frame is dr dt F =dp dt +CI Fde dt F +CI F(ωF×e)(4.25) 4.3 Moving Frame Analysis 73 If (4.25) is multiplied by the rotation matrix defined in (4.20), (4.25) becomes CF Idr dt I =dp dt F +de dt (4.26) The velocity components at point R are dr dt = [˙x˙y˙z]T . The error components are defined as e= [eweleo]T , being the forward, lateral and longitudinal errors respectively, and the term ωF×ecan be expanded to ωF×e=   −˙ χfsinγf ˙ γf ˙ χfcosγf   ×   ew el eo   =   ˙ γfeo−˙ χfcosγfel ˙ χfcosγfew+˙ χfsinγfeo −˙ χfsinγfel−˙ γfew   (4.27) . From (4.26), and using the previous definitions, it can be rewritten as CF I   ˙x ˙y ˙z   =   ˙s+˙ew+˙ γfeo−˙ χfcosγfel ˙el+˙ χfcosγfew+˙ χfsinγfeo ˙eo−˙ χfsinγfel−˙ γfew   (4.28) Solving for the error rate of change components yields ˙ew=Vgcosγ−γfcosχ−χf−˙s˙ γfeo+˙ χfcosγfel(4.29a) ˙el=Vgcosγ−γfsinχ−χf−˙ χfcosγfew−˙ χfsinγfeo(4.29b) ˙eo=Vgcosγf−γ+˙ χfsinγfel+˙ γfew(4.29c) Note that the course and groundspeed are used because ground-reference measurements produces a better response in the presence of stochastic wind disturbances, such as gusts [72]. The COTS autopilot can implement course and flight path-hold loops, the course and flight path angle dynamics can be expressed as: ˙ χ=ς(χc−χ)(4.30a) ˙ γ=ϑ(γc−γ)(4.30b) where γc and χc are the commanded flight path and course angles and ς and ϑ are constants that characterize the speed of response of such angles. The angular error is 74 Chapter 4. Wind efficient flight and trajectory tracking defined as: ˙ χe=˙ χf−ς(χc−χ)(4.31a) ˙ γe=˙ γf−ϑ(γc−γ)(4.31b) The tracking problem can be formulated from (4.30) and (4.31) , in which the goal is to minimize each one of these equations until the values reach a point close to zero. Considering the development on the vector field based curve tracking presented by Nelson and Miao [ 72 , 48 ], the flight path angle errors can be expressed as a function of the positional error (4.30): χl(el) = −χ∞ 2 πarctan(κel)(4.32a) γo(eo) = −γ∞ 2 πarctan(µeo)(4.32b) where −χ∞∈0,π 2 and −γ∞∈0,π 2 are the maximum absolute desired value of the respective error. From here, the Lyapunov function L is proposed with three objectives, one is to minimize the distance error as close as possible to zero, then is to minimize in the same way both the course and the flight path errors. L=1 2νe2 w+e2 l+e2 o+1 2(χ−χl(el))2+1 2(γ−γo(eo))2(4.33) where ν is a weighting factor, that reduces the contribution of the positional error. The three terms are defined for the minimization of the different errors with the aforementioned dynamics. The time derivative of the Lyapunov function in (4.33) is: ˙ L=1 ν(˙ewew+˙elel+˙eoeo)+(χ−χl(el))( ˙ χ−˙ χl(el))+(γ−γl(eo))( ˙ γ−˙ γl(eo)) (4.34) 4.3 Moving Frame Analysis 75 Substituting (4.29) into (4.34) yields: ˙ L=1 νVgcosγ−γfcosχ−χf−˙s˙ γfeo+˙ χfcosγfelew +Vgcosγ−γfsinχ−χf−˙ χfcosγfew−˙ χfsinγfeoel +Vgcosγf−γ+˙ χfsinγfel+˙ γfeweo +(χ−χl(el))˙ χf−ς(χc−χ)−˙ χl(el) +(γ−γo(eo))˙ γf−ϑ(γc−γ)−˙ γo(eo) (4.35) Rearranging the equation (4.35) and solving for γc and χc , the control laws for flight path angle and track are determined: χc(el) =Kχ(χ−χl(el))+ χ+1 αχl(el)Vgcosγsin χ −1 α elVgcosγ(sin χ−sin χl(el)) ν(χ−χl(el)) (4.36) γc=−Kγ(γ−γl(eo))+γ+1 β˙ γo(eo)Vgsinγ −1 β eoVg(sinγ−sinγo(eo)) ν(γ−γo(eo)) (4.37) where Kχ , Kγ and ν , α and β are the coefficient that make the state variables tend to zero. To proof the stability of the control law, substituting the (4.36) and (4.37) into (4.34) , the equation can be rewritten as: ˙ L=1 νelVgcosγsin(χl(el))−1 νeoVgsin(γo(eo))−Kχ(χ−χl(el))2−Kγ(γ−γl(eo))2 (4.38) It is easily observed that when the errors tend to zero, χd=χf and γo=γf , then ˙ L=0 and if any of these conditions is not truth then ˙ L<0 . When time tends to infinity, all the errors, including the values of γo(eo)and χl(el), tend to zero. The most important advantage of this adaptive control method is that it is completely independent of the dynamic model of the UAS, which is very helpful for uncommon platform configurations, such as flying wing, or delta wing UAS’s. In addition, the error equations are relatively easy to compute and the implementation will not affect the overall 76 Chapter 4. Wind efficient flight and trajectory tracking -3000 -2500 -2000 -1500 -1000 -500 0 500 1000 -1000 -800 -600 -400 -200 0 200 400 600 800 1000 350 400 450 500 x (m) y (m) height (m) flight trajectory with x-y-height referened path with x-y-height start No.13 No.12 No.8 No.9 No.11 No.10 No.7 No.5 No.6 No.1 No.2 No.3 No.4 No.14 Figure 4.9 As presented in [ 48 ], the experiment shows convergence in the path marked in black for a total of 14 waypoints and an offtrack starting position. computational time of the other functionalities presented. In [ 48 ] the algorithm was able to follow accurately a 3D path (see Figure 4.9), however they claim that the system is able to follow a 5th order spline without any experimental proof of it. It is assumed that the main reason of the selected flight path of experiment is the easy implementation of the course and flight path angle of the desired trajectory for the autopilot control. In Chapter 6, the software of the selected autopilot for this research is analyzed since the actual control loop of the groundspeed and course may require additional implementations. The adaptive method has been proved to work in straight trajectories as in [ 72 ], the stability conditions for both the error dynamics and the control laws are analyzed extensively. Furthermore, in [ 48 ], the method is proved at a experimental stage for straight 3D paths with great convergence above. There has been attempts to prove the vector field theory approach in more complex situations at an experimental phase. Nelson [ 72 ] simplifies the problem by approximate the turns to arcs, however Jung [ 74 ] considers a more complex curve for the 2-Dimensional approach. However the results achieved show that with low wind disturbances and with a considerable level of smoothness along the path, the control shall be able to track the splines effectively. 4.4 Chapter summary and conclusions 77 4.4 Chapter summary and conclusions This chapter discussed two processes. The sequence reconfiguration and the trajectory generation, which are the primary methods for the realization of the flight duration increase. In the case of the sequence reconfiguration, an algorithm was proposed. This algorithm reconfigures the waypoint sequence, which was previously generated with an area partition method, by modifying the weight of the cells of a partitioned area according to the wind information; this information can be provided to the system either with the use of an external sensor, such as the anemometer, or while flying to perform an online reconfiguration based on the determination of the predominant wind speed. This reconfiguration intends to maximize the legs in which the aircraft flies with a tailwind. This provokes an automatic reduction of the thrust and an energy saving. The sequencing reconfiguration was incorporated late into the design of the system for flight duration increase. Initially, the goal was to create a full flight plan that uses the trajectory generation algorithm, so that the tracking is performed at every leg of the flight plan. After an analysis of the complexity the tracking algorithm it was decided that this function could be triggered only at specific locations of the flight plan, which are previously analyzed, and reduces the complexity of the tracking functionality. The trajectory tracking process focuses on the ability of the system to follow the Rayleigh cycles, which is accomplished after a previous selection of legs that are candidates for dynamic soaring. First, a moving frame is added to the curve, based on the Type-1 Bishop frame, that determines the attitude that the aircraft has to follow while flying in a tracking point along the curve. Then, after the framing process, an adaptive control system, based on the vector field theory is implemented. An error calculation is performed in which the positional and angle errors are determined from the current aircraft position with respect the frame at a given tracking point. This error leads to the construction of a Lyapunov function, that determines the control laws that will minimize this function and ultimately follow the desired trajectory. The main contribution discussed in this chapter is the inclusion of a waypoint sequencing algorithm to loose up the complexity of tracking continuous Rayleigh cycles that may produce positional errors on a given mission. This allows the utilization of a two-stage control strategy that uses the legacy autopilot control for waypoint navigation, which typically uses a classic outer guidance loop and an inner control loop, and also the implementation of the vector-field adaptive control discussed above. The next chapter discusses the integration of the wind-based waypoint sequencing, and the Rayleigh cycle tracking algorithms with the other functionalities that were discussed 78 Chapter 4. Wind efficient flight and trajectory tracking in the previous chapters. Also, the different testing benches are discussed. 5 System integration Foolproof systems don’t take into account the ingenuity of fools. Gene Brown Th e previous chapters explained the different functions that are necessary to extend the flight duration for oceanic and coastal missions. The main functionalities that are considered for the duration extension are: Wind field identification function : identifies the wind velocity vector and the wind velocity rate of change at the given location. Wind feature characterization function : based on the available wind information, i.e. wind estimates and weather information, some features of interest are identified and characterized to determine if the wind field is suitable for AEH. Trajectory generation function : aims to generate suitable trajectories on a waypoint- to-waypoint basis to perform Rayleigh cycles. Waypoint sequence reconfiguration function : given a waypoint sequence, some legs are selected in order to perform the Rayleigh cycles, and the sequence is reconfigured for maximum energy gain. Trajectory tracking function : in the legs that are selected for Rayleigh cycles, the generated trajectory has to be tracked with a robust control scheme for a feasible energy gain. Based on these functions, an architecture needs to be defined to implement them into a UAS. This architecture has to fulfill a set of requirements for robustness, safety and reliability to be able to provide with real-time flight duration enhancement. This architecture shall meet stringent computational cost, due to the complexity and has to allow upgradeability and maintainability. With this philosophy, this chapter proposes an architecture of a Flight Duration Enhancement System (FDES) byfollowing a design methodology of systems engineering that will ensure the correct functionality of the system. 79 86 Chapter 5. System integration Wind Estimation Wind Vector Wind Dir. Wind Accel. Feature Identification Shear params. Gust params. Wind map Trajectory Generation Wpt seq. Rayleigh legs Trajectory Tracking Curve frames Error tracking Control cmds. Communication framework (Topic Subscribe / Publish Scheme) External Communications Interface (Sockets - ROS/MAVLINK) Autopilot Figure 5.4 The Software Architecture defined, focused on the communication framework scheme. The low level functions publish and subscribe to topics to perform their different routines in an asynchronous manner. processes the data if there are subscribers connected on the created topic. A simple example of the publishing is presented in Code 5.1. Code 5.1 Simple publishing example. ros::NodeHandle n; ros::Publisher pub = n.advertise<std_msgs::String>("name", 5); std_msgs::StringPtr str(new std_msgs::String); pub.publish(str); The subscriber method is easier to implement, as the routine subscribe() , of the needs to be called either by a a class, a function or a functor object, which leads to a change in the way the method is called. Code 5.2 Simple subscribing function example. void callback(const std_msgs::StringConstPtr& str) { 5.2 Software integration 87 ... } ... ros::Subscriber sub = nh.subscribe("my_topic", 1, callback); At each function of the FDES, there are a series of publishers and subscribers that transfer the data from each function to the other. The processes are executed asynchronously, and the trigger of the execution of the functionalities is the existence of the different topics that the routine is subscribing to, e.g. if the wind calculation requires information from the sensors, it waits if the message are available in MAVLINK and performs the calculation and publishes the necessary topic. Once the communication framework is implemented, the next part is the definition on how the different functionalities were implemented. 5.2.2 Wind field estimation and feature identification As discussed in Section 2.6, the wind field estimation (wind velocity and wind velocity rate of change) depends on different values that reside in the autopilot. The necessary messages, for the wind field estimation and the feature identification algorithm that the autopilot transmits to the FDES trough the communication framework are enumerated below: 1. Gyro status defined in the MAVLINK message MAV_SYS_STATUS_SENSOR_3D_- GYRO message of the MAV_SYS_STATUS_SENSOR type enumeration. 2. Acceleration status defined in the MAVLINK message MAV_SYS_STATUS_SENSOR_- 3D_ACCEL of the MAV_SYS_STATUS_SENSOR type enumeration. 3. Accelerations defined in the xacc , yacc , zacc MAVLINK messages of the SCALED IMU type enumeration. 4. Rotational rates defined in the xgyro , ygyro , zgyro MAVLINK messages of the SCALED IMU type enumeration. 5. Attitude angles defined in the roll , pitch , yaw MAVLINK messages of the ATTITUDE type enumeration 6. Groundspeed , defined in the vx , vy , and vz messages of the GLOBAL_POSITION_- INT type enumeration. 88 Chapter 5. System integration 7. Covariance matrix , defined in the covariance message of the GLOBAL_POSITION_- INT_COV , additionally, the velocity covariance in the velvariance message of the CONTROL_SYSTEM_STATE type enumeration. 8. Airspeed , defined in the airspeed message of the CONTROL_SYSTEM_STATE type enumeration. Once the wind parameters are identified, the computational functions of the wind functionality of the FDES are implemented following the algorithms of Section 2.6. The result of these functions rewrite the messages of the autopilot that are related to the default wind calculation and publishes new messages in the MAVLINK format related to the identification of the wind. 5.2.3 Trajectory generation and sequence reconfiguration implementation The trajectory generation and sequence reconfiguration function resides in the FDES module. The implementation requires of a valid waypoint lists and valid wind information. As it is mentioned before, the minimum information necessary to implement this functionality is the horizontal predominant wind direction, this information may be extracted from meteorological reports, if an Internet connection is available. Also, this information can be directly input by the user, or can be obtained with actual measurements of the UAS, since the running mean and standard deviation is stored in the system. If the information is used pre-flight, the system performs the necessary calculations and proposes a new sequence of waypoints from the beginning, and if the information is onboard, then the system proposes a reconfiguration of the remaining waypoints. In addition, the system generates candidate legs for the implementation of Rayleigh cycles, and sends that leg information to the trajectory generation function that is performed in two stages, one that determines the auxiliary waypoints of that leg, i.e. the control points of the curve, including the generation of auxiliary waypoints to implement a full cycle or a portion of it. The following diagram illustrates the process described above. 5.2.3.1 Curve framing implementation The curve framing was implemented as a separated module that receives the next Rayleigh leg from the current waypoint. The function performs the framing using the Type 1 Bishop frame, firstly by determining the number of tracking points among the curve which depend on the distance, typically a tracking point will be at approximately every meter of the leg. Once the curve is generated, the tracking points and the unitary vector associated to the frame is stored in a vector type variable. Once the next waypoint 5.2 Software integration 89 Wpt Sequence Wind information: * Wind direction * Wind magnitude * Shear parameters Sequence Reconfiguration Determination of Rayleigh Legs Determination of auxiliary points Trajectory calculation function New Wpt. Sequence Coordinates of legs Calculation of velocity fcts. Figure 5.5 Functional diagram of the trajectory generation function. is the beginning of a Rayleigh leg ( MAV_ROI_WPNEXT of the MAV_ROI type enumeration), the curve framing implementation triggers the tracking function that is described in the following section. 5.2.4 Control Implementation The control implementation is implemented as a selection function that will enable the tracking of the mission with the autopilot algorithms, which determine the position error and generate the commands using a PD controller. At this stage, the addition of this research is the tuning parameter selection, that is modified based on the wind. Observations on the extensive test campaign helped to determine that the tuning can be changed depending on the wind phenomena, in the case of heavy winds, shear or gusts, the system will lower the gains, and in the presence of calmed winds the system will tighten the gains to ensure less error. The preset values of the tuning are stored as constant values, and the system will choose a preset value based on the environmental conditions. As it is mentioned in the previous section, if there is a candidate leg to execute a Rayleigh cycle, the auxiliary waypoints and the control points were generated and the curve is framed, then the system will start calculating the position error between the first tracking point. This is the beginning of the Rayleigh cycle. Once this error is minimum, which will depend on the acceptance criteria that is modified based on the environmental conditions, then the control laws will be determined and will be translated into pitch, roll and yaw commands to track the curve. At every step, the system is monitoring the distance to the tracking waypoint and it will determine, based on the heading and the position, if the tracking waypoint is reachable or not. If it is not reachable, the function will switch to the next waypoint, and the process will start again. 90 Chapter 5. System integration To prevent latency with the communications between the FDES the trajectory error determination and the tracking function were implemented as functions of the PX4 firmware, meaning the development framework of the PX4 was utilized [ 75 ]. The implemented function receives the tracking waypoint and the Bishop frame unitary vectors, computes the error equations presented in (4.30) and (4.31) , and determines the control commands while monitoring for a change of the tracking waypoint which is computed in the FDES module. 5.2.5 Flying modes and control parameters The PX4 architecture contains different flight modes that provide the user with options on the interaction level between the pilot and the FDES. To fully implement the FDES in a safely manner, two additional modules were added to the preexisting modes. The two added modes are: RAYLEIGH . This mode enables the reconfiguration of the flight path and the selection of legs for dynamic soaring. The tracking function is fully enabled allowing the dynamic change of the control gains. WIND . This mode enables the reconfiguration of the flight path only. The tracking function is disabled, but the dynamic change of the control parameters is allowed. The control parameters that have to be modified for tuning are the feed forward gain, the rate gain and the trim offsets with the integrator gain for the roll, the pitch and optionally the yaw. Furthermore, the time constant of the outer loop, or guidance loop is also adjusted for both the pitch and the roll response. Regarding the waypoint tracking, which is not part of the Rayleigh cycle tracking, the controller is based on an L1 adaptive control technique [ 77 ], which has different tuning parameters that can be changed. The most critical parameters are the L1 damping factor the control and the L1 period, which is the distance that defines the tracking point ahead of the aircraft that is following. Finally, the last set of parameters that have to be tuned are those related to the Total Energy Control System (TECS), that aim to control the vertical response of the aircraft when it is climbing, descending or in level flight. The details of how the parameter values for high and low winds were determined can be found in Section 6.3.3. 5.3 Testing benches The validation and verification is a key objective in the design of the FDES. Therefore, some testing benches need to be carefully designed to validate the different functionalities 5.3 Testing benches 91 at various levels. For this purpose, two main approaches were selected. The first one consists in the creation of a Hardware-In-The-Loop (HITL) architecture, that aims to prove the communications and the reliability of the hardware. It allows the testing of the different flight modes before the actual field experiments, but also to test the behavior of the systems under some potential failures, which is a significant step in the validation phase. The Software-In-The-Loop (SITL) architecture aims to test the different functions of the FDES at various levels and with different methodologies. It is key for testing modifications on the functionalities that are necessary to implement the control routines in the autopilot and the other functions in the FDES module. The software shall allow the testing with a dynamic model but also with available data from other flights, such as the estimation of the wind field, the identification of features, the control responses, etc. This section describes the architecture of both testing benches and discusses the selection of such architectures and the impact on the system validation. 5.3.1 Hardware-In-The-Loop setup The HITL bench is composed by three main hardware components: The autopilot, that has updated with the latest version of the firmware, and has the telemetry module and the RC system attached to it. The FDES module, that has a USB-to-UART adapter for autopilot communications, has the ROS. The ground-station computer that runs a flying simulation software and the groundstation software that receives the data from the telemetry module. For the flying simulation software, a commercial license of XPlane-10 was utilized as it is the main supported platform for the PX4 developer team. The communication is achieved with UDP ports that are utilized for the XPlane to transmit the flight parameters. The Simulator runs the dynamic model of the aircraft, that is achieved by compiling the JSBSim library [ 78 ], which is is an open source flight dynamics model (FDM). The simulator stablishes direct communication with the ground station, that runs the QGround- Control software [ 79 ] and embeds the simulator parameters (position, velocity, etc.) into MAVLINK messages with a unique message type enumeration. The Qgroudncontrol is connected trough a USB-To-UART interface to the autopilot that computes the mission parameters and sends the commands that will be embedded into MAVLINK command messages. The command parameters are sent to the simulator which provides a real-time response and offers the versatility of testing different flight modes. Since the autopilot has 92 Chapter 5. System integration Figure 5.6 The HITL functional architecture, in which the ground station is connected to the flight controller hardware via USB, and acts as a gateway to forward data between the X-Plane simulator running on a development computer, PX4, and any offboard API. the telemetry module and the RC system attached to it, the bench allows pilot commands from the RC transmitter and also the possibility of getting a whole set of telemetry parameters into another ground station for communication testing. The HITL architecture diagram is shown in Figure 5.6. 5.3.2 Software-In-The-Loop and simulation bench setup The SITL bench is implemented in two parts. The first part consist on the generation and testing of the different models and functionalities that are utilized in a simple simulation form. This is accomplished with MATLAB ® SIMULINK ® which allows the creation of complex atmospheric scenarios that will allow the testing of the different functionalities of the FDES before the full implementation on the system. The Simulink Aerospace Toolbox allows the implementation of shear, gusts, and different atmospheric conditions, such as the humidity, dew point, temperature, etc. Once the complex scenario is created, the system allows the simulation using MATLAB ® S-Functions, which enable the import of the libraries that were created for the FDES. The second part aims to test different changes on the source code of the PX4 autopilot. As it is mentioned before, the trajectory generation function, the flight modes and the created MAVLINK messages require the modification of the source code and the creation of a firmware image particular for the implementation of the system. Architecture-wise the SITL is similar to the HITL architecture illustrated in Figure 5.6, however it does not require any hardware and can be fully implemented into one station. The purpose of the 5.4 Chapter summary and conclusions 93 Figure 5.7 The SITL bench architecture shows the different modules and the communication means of the bench. SITL test bench allow rapid prototyping and testing of different functionalities before compiling it into a final version. Both the HITL and SITL test benches allow a full functional testing of the aforementioned architecture before testing the system on-board. The off-board testing was a very important part of the design of the FDES, which permits that different configurations and software versions at various levels are being explored. 5.4 Chapter summary and conclusions This chapter discussed the integration of the different functionalities into a Flight Duration Enhancement System (FDES), that follows a systems engineering approach to implement the software and hardware architecture. A set of system requirements are defined, from which the system architecture is derived, and then cascaded into the implementation of the different functionalities. In this chapter, all the design decisions, such as the selection of hardware components, the definition of the software and hardware architectures, etc. are justified with a based on the defined requirements. The design philosophy aimed to produce a system that can be scalable, and has high levels of maintainability and upgradeability without modifying the components or other pieces of software. In terms of the hardware architecture the main design decision was the utilization of a FDES module with powerful characteristics in order to implement the different algorithms that were discussed in the previous chapters. Also, the selection of a COTS autopilot and sensors allow the implementation of the different functionalities as expected. 94 Chapter 5. System integration The software architecture had different design decisions that include the selection of a loosely coupled architecture in which the different functionalities are enabled and communicated with a topic-publish and subscribe scheme. In addition the selection of the means for interaction between the functionalities, the selection of the main communication protocol were discussed and justified according the system requirements. Finally, the definition of the Software-In-The-Loop (SITL) and Hardware-In-The-Loop (HITL) architectures were defined. These functions are the foundation of the testing activities that will determine if the different hypothesis presented throughout this thesis are correct. The main contribution of this chapter was the utilization of a systems-engineering approach in order to define and implement in a well-defined structure the different functions. The systems-engineering model is not a typical approach in research, however, its utilization permits the definition of appropriate test to validate and verify the different functionalities at various levels. Therefore, each functionality could be tested separately even at a prototype phase, and also during integration, allowing the utilization of real-telemetry data, or the testing of new missions. The next chapter describes the different tests that were performed to each functionality at the stages described above. Conjointly, the description of the different testing platforms and sites is presented showing the results of the extensive experimental campaign that took place in order to validate the system. 6 Systems validation and verification, real experiments and results analysis It works! It works! Dr. Emmet Brown Th is chapter describes the experimental process followed throughout the design and implementation of the FDES. Since the design approach followed the systems engineering methodology, there was an extensive testing campaign throughout the design and implementation process. These testing stages intend to validate each function at high an low level. The testing process started at the prototyping of the implementation of the different functionalities, and for this the MATLAB ® SIMULINK ® environment was utilized. Once the functions were properly tested in the simulation environment, the next stage was to implement it in the actual hardware followed by the SITL and HITL testing. After this stage, the system was tested in different locations and scenarios including the metropolitan Seville region, the coastal region of the South of Spain and some experiments that were performed in collaboration with other partners of the MarineUAS project. The chapter shows the validation and verification methods employed at each stage of the development, including a description of the experimental platforms and the insights of the different testing methods. 6.1 Wind field estimation and features identification testing The first part of the development was the creation of algorithms for the estimation of the wind field (wind velocity and wind velocity rate of change). For this purpose, the initial stage was the implementation of the wind models in MATLAB ® SIMULINK ® environment, together with the implementation of the dynamic model presented in Chapter 2, with some 95 102 Chapter 6. Systems validation and verification, real experiments and results analysis 0 50 100 150 200 250 300 350 400 Time (s) 0 0.5 1 1.5 2 2.5 3 Wind Speed (m/s) (a) 0 50 100 150 200 250 300 350 400 Time (s) -200 -150 -100 -50 0 50 100 150 200 Wind Direction (deg) (b) Figure 6.7 The results obtained with the autopilot calculation of the wind speed and direction on the experiment presented in Figure 6.6a, in which (a) shows the wind speed and (b) show the wind direction, in which abrupt changes are observed despite the assessment on the region stability and roughness. The Weibull probability plot that illustrates the dataset can be seen in Figure 6.10, proving that the dataset matches a Weibull distribution, thus, the results are statistically coherent. As it is mentioned before the average wind speed is in the range of the wind speed reported in the Airport station. Taking the whole dataset, the algorithm intended to find a characterization curve of the wind-shear. Initially the experiments intended to fit a simpler model of the shear, based on an empirical power law, which is described in (6.1). W1 W2 =h1 h2 ε (6.1) The goal was to find a coefficient ε , namely the Prandtl exponent, that provides with 6.1 Wind field estimation and features identification testing 103 101102 Data 0.001 0.003 0.01 0.02 0.05 0.10 0.25 0.50 0.75 0.90 0.96 0.99 0.999 Probability Figure 6.8 Weibull distribution plot of wind speed magnitude calculated by the PX4 in the helical portion of the flight shown in Figure 6.6a. TIme (sec) 0200 400 Wind Speed (m/s) -25 -20 -15 -10 -5 0 5 10 15 20 25 (a) Time (sec) 0200 400 Wind Speed (m/s) -10 -5 0 5 10 15 (b) Time (sec) 0200 400 Wind Speed (m/s) -30 -20 -10 0 10 20 30 (c) Time (sec) 0200 400 Wind Speed Error 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 (d) Figure 6.9 The estimation results with the direct computation method: (a) shows the north component, (b) corresponds to the east component, (c) is the down component, and (d) is the estimated RMSE error on the total magnitude. the best possible fit of the dataset. The model evolved in later experiments, allowing the characterization of a more precise model, and simplifying the fitting process, since the goal became to find the coefficient Aσin (2.9), which becomes into a simple problem. 104 Chapter 6. Systems validation and verification, real experiments and results analysis 100 10 1 Data 0t001 0t00. 0t01 0t03 0t05 0t10 0t35 0t50 0t75 0t90 0t96 0t99 0t999 Probability Figure 6.10 Weibull distribution plot of wind speed magnitude calculated off-board with the direct computation method in the helical portion of the flight shown in Figure 6.6a. The red dots indicate the theoretical Weibull distribution with the proposed parameters and the blue dots represent the dataset of the function. 810 12 14 16 18 20 22 24 26 28 Wind (m/s) 0 20 40 60 80 100 120 140 160 Altitude (m) Figure 6.11 Determination of Prandtl exponent with the data estimates obtained with all the dataset of the flight Figure 6.6a. The results obtained with the empirical power law can be observed in Figure 6.11. It is suggested that the value of ε=0.14 is applicable to a wind shear, however it was observed that a value of ε=0.18 was more accurate for this particular dataset. The prediction process evolved from the determination to the Prandtl coefficient ε to the determination of the aforementioned Aσ coefficient, which allows a more computational efficient process. The Prandtl coefficient requires different comparisons among wind speeds to fit in the best possible curve, however, in the case of the estimation of Aσ , the process is simplified as only two reference values are required. The results of the glider in the flight presented in Figure 6.6b followed the same estimation process, however, the use of the empirical power law was substituted to the use of equation (2.9) . The data was fitted to a curve, in which the resultant coefficient is 6.1 Wind field estimation and features identification testing 105 -20000 -15000 -10000 -5000 0 5000 X Position (m) -8000 -6000 -4000 -2000 0 2000 4000 6000 Y position (m) Scenario 1 and 2 Scenario 3 and 4 Scenario 5 and 6 Figure 6.12 Trajectories of the simulated flights that show the effect of the different wind features when the aircraft is subject to them and the controller is unable to reject the disturbance and loses the intended trajectory track. In these simulations, the heading control was enabled and the waypoint control was disabled. Aσ=1.43. With the evolution of the prediction process, the function was re-tested in a simulated scenario, which consisted on testing different cases that combine the wind features. The tested flight on a squared trajectory, with the implementation of the waypoint tracking algorithm of the PX4. The flights have 2 types of changes in altitude, one that has a sinusoidal wind vertical profile and the other a stepped profile with random altitude changes and holds. The MATLAB ® Simulink ® environment was utilized for this simulation and the simulation characteristics are shown in Table 6.4. Table 6.4 Simulation 2 Scenarios Characteristics. Scenario Wind Features ∆Altitude 1 Constant Wind High (Sinusoidal) 2 Constant Wind Low (Steps) 3 Wind Shear High (Sinusoidal) 4 Wind Shear Low (Steps) 5 Wind Shear + Gust High (Sinusoidal) 6 Wind Shear + Gust Low (Steps) Each scenario was simulated for 1500 seconds, and a maximum shear speed was considered to be Wshearmax =15 at h=300 m. The tested horizontal profiles are shown in Figure 6.12 and the horizontal profiles are depicted in Figure 6.13. 106 Chapter 6. Systems validation and verification, real experiments and results analysis 0 500 1000 1500 Time (s) 0 50 100 150 200 250 Altitude (m) Scenarios 2,4,6 Scenarios 1,3,5 Figure 6.13 Vertical profiles of the simulated flights of simulation 2, showing the stepped vertical profile in blue and the sinusoidal profile in red. The altitude changes allow the estimation of the wind field in a broader altitude spectrum, which enhances the prediction process. The simulation intends to validate the ability of the system to react to different features and to analyze the convergence time of the prediction. At every time, the prediction process intends to estimate the best fit to characterize the shear layer. The goal of the first two scenarios is to present very low winds to detect if the system rejects the noise and does not consider it a feature, even with the presence of an altitude change. The third and fourth scenarios have only a shear model activated with a different vertical profile. The last two cases consider the presence of a discrete gust that appears after a few minutes after the UAS is airborne. The results of the estimation include the average wind speed µ(W) , the wind-shear coefficient A , and the magnitude of a discrete gust if it is present Wm as in Section 2.1.1.1, which is sufficient to assess the efficiency of the wind-estimation and identification functions. Table 6.5 Scenarios outputs (single run). Scenario µ(W)A Wm 10.1695m/s N/A N/A 20.1694m/s N/A N/A 311.2979m/s 1.5384 N/A 411.2972m/s 1.7982 N/A 515.5237m/s 1.4415,1.8129 4.96m/s 615.5176m/s 1.2712,1.1321 4.94m/s 6.1 Wind field estimation and features identification testing 107 6 7 8 9 10 11 12 13 14 Wind Speed(m/s) 0 50 100 150 200 250 300 Altitude (m) Instantaneous wind estimations Predicted wind profile at 30 seconds Predicted wind profile at 50 seconds Predicted wind profile at 60 seconds Predicted wind profile at 70 seconds (a) 6 8 10 12 14 16 18 20 0 50 100 150 200 250 300 Instantaneous wind estimations Predicted wind profile at 30 seconds Predicted wind profile at 70 seconds Predicted wind profile 30 seconds after gust Predicted wind profile 50 seconds after gust Predicted wind profile 70 seconds after gust Predicted wind profile 100 seconds after gust Wind Speed(m/s) Altitude (m) (b) 5 6 7 8 9 10 11 12 13 14 0 50 100 150 200 250 300 Instantaneous wind estimations Predicted wind profile at 30 seconds Predicted wind profile at 40 seconds Predicted wind profile at 50 seconds Predicted wind profile at 60 seconds Wind Speed(m/s) Altitude (m) (c) Figure 6.14 The evolution of the prediction process presented in the scenarios shown in Table 6.5, considering only the scenarios with the presence of a wind shear on sinusoidal and stepped vertical profiles. (a) and (c) correspond to the scenario 3 and 4 respectively, in which it can be observed that the convergence time is faster in a more varying vertical profile. (b) shows the scenario 5, in which a discrete gust provokes a sudden change on the estimation region. As it is mentioned before, the prediction process is continuously evolving with the addition of new data into the data-set, and even the inclusion of external data allow the possibility of identifying the features more accurately. This process evolution in time is illustrated in Figure 6.14 in which the development of the prediction towards the final convergence value is illustrated. The goal was to determine the convergence time, to calculate the computational time of the identification cycle to assess the impact on the overall computational time and also to prove the ability to detect multiple features. The results presented in Figure 6.14 proved that the method of looking for jumps in the standard deviation is adequate to detect external features, such as, discrete gusts. In the case of the computational cost, the evolution process has a GA algorithm which may impact significantly the overall process of flight duration enhancement. However, by utilizing the clocks function of c++ time.h library, the computational time can be determined by dividing the result into the constant CLOCKS_PER_SECOND which is part of the same 108 Chapter 6. Systems validation and verification, real experiments and results analysis Table 6.6 SkyWalker Characteristics. Wing Span 2122 mm Wing Area 80 dm2 Max Payload 2 kg Center of Gravity 435 mm away from nose (a) (b) Figure 6.15 The SkyWalker X8 (a), from (SkyWalker Technology Co. Ltd, Wuhan, China) is a flying wing platform with a carbon fiber frame. The 4525DO sensor (TE Connectivity Ltd., Shaffhausen, Switzerland) with a resolution of 0.84 Pa [80]. library. The result shows an average execution time of µ(tGA) = 0.053s with an standard deviation of σ(tGA) = 0.012s. 6.1.2 Wind estimation and feature identification SITL tests After the wind field estimation and wind field identification function were validated at a prototype level, the next stage was to prove the estimation and identification process in a SITL scenario, this occurred after the creation of header files and implementation file. The implementation uses a Simulink ® that interacts with the ground station software to replicate entirely a series of flights that were performed in the region of Brenes, Seville with another platform, the SkyWalker X8 (see Figure 6.15), which is a very popular commercial platform due to its payload capabilities. Its characteristics are enumerated in Table 6.6. The relevant information of the experiments that were replicated in the SITL test can be observed in Table 6.7, and the trajectories of the corresponding flights can be observed in Figure 6.16. 6.1 Wind field estimation and features identification testing 109 37.524 37.523 37.522 37.521 Latitude (deg) 37.52 37.519 37.518 37.517 -5.862 -5.86 Longitude (deg) -5.858 -5.856 40 20 60 80 100 120 140 160 180 -5.854 Altitude (m) (a) 37.522 37.521 37.52 Latitude (deg) 37.519 37.518 37.517 -5.86 -5.859 -5.858 Longitude (deg) -5.857 -5.856 160 20 40 60 80 100 120 140 -5.855 Altitude (m) (b) 37.524 37.523 37.522 37.521 Latitude (deg) 37.52 37.519 37.518 37.517 -5.861 -5.86 -5.859 Longitude (deg) -5.858 -5.857 -5.856 160 20 40 60 80 100 120 140 -5.855 Altitude (m) (c) 37.523 37.522 37.521 Latitude (deg) 37.52 37.519 37.518 37.517 -5.861 -5.86 -5.859 Longitude (deg) -5.858 -5.857 -5.856 160 20 40 60 80 100 120 140 -5.855 Altitude (m) (d) Figure 6.16 The trajectories of the flights performed by the SkyWalker X8 in the region of Brenes, Seville. The trajectories show helical flights with at altitude that goes from 100m to 150m above ground. Table 6.7 Experiments information. Flight 1 Flight 2 Flight 3 Flight 4 Duration 521 s 315 s 631 s 749 s Distance Traveled 5.1 km 3.7 km 6.3 km 7.4 km Maximum Altitude 179 m 125 m 134 m 146 m The estimation results obtained from the four flights have a more noisy input, which is attributed with the windy days with low stability in the region at the time of the flights. Nevertheless, the system was able to perform an estimate with these noisy measurements, excluding the outlying points that were caused due to the instability of the system. The wind vs. altitude maps of the flights can be observe in Figure 6.17 Note that in Figure 6.17b the results show two curves, which is similar to the simulation result showed in Figure 6.14b in which a jump in the standard deviation was interpreted as a discrete gusts, and then the curve is characterized again. The scenarios allow understanding how the noise in the estimation, derived from the airspeed measurement and the region stability affect the prediction process and the estimation. In the wind estimation data 110 Chapter 6. Systems validation and verification, real experiments and results analysis Wind Speed (m/s) 0 2 4 6 8 10 12 Altitude (m) 20 40 60 80 100 120 140 160 180 Wind Estimation Most Probable WS Wind Prediction (last) (a) Wind Speed (m/s) 0123456789 ALTITUDE (m) 40 60 80 100 120 140 160 Wind estimations Most probable WS Prediction (b. gust) Prediction (a. gust) (b) Wind Speed (m/s) 0 2 4 6 8 10 12 Altitude (m) 20 40 60 80 100 120 140 160 Wind Estimations Most Probable WS Wind Prediction (c) Wind Speed (m/s) 0 2 4 6 8 10 12 14 Altitude (m) 20 40 60 80 100 120 140 160 Wind Estimation Most Probable Wind Speed Wind Prediction (d) Figure 6.17 The estimation, most probable wind speeds and last shear prediction of the four flights presented in Figure 6.16. The blue points represent the wind speed estimations performed with the direct computation method, the orange dots represent the evolution of the most probable wind speed at the respective altitude and the lines represent the latest shear curve. observed in Figure 6.17b, Figure 6.17c, and Figure 6.17d it can be seen a concentration of disperse wind estimations in the altitude in which the helical trajectory is performed. The concentrated data is distributed following the Weibull model allowing at a more accurate data prediction at those stages. Moreover, the actual variability of the data may not present a significant impact on the duration enhancement process, since the goal is to find just adequate conditions for the execution of the soaring maneuver and not to get a precise set of measurements. 6.2 Trajectory-related simulations In the case of the trajectory generation and sequence reconfiguration, there were different experiments that were considered in to validate these functionalities. These experiments consider the use of the SITL testing bench described in the previous chapter, including also the software developed by Balampanis [ 67 ] for partitioning complex areas and the generation of an initial waypoint sequence, and also the use of telemetry data for afterprocessing to test the different capabilities. 6.2 Trajectory-related simulations 111 Figure 6.18 Test case area in the north of the Seattle Metropolitan Area in the state of Washington in the United States. The red-line polygon defines the mission area. X position (m) #104 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 Y Position (m) #104 -5 -4 -3 -2 -1 0 1 2 Starting Point DIrection of sequence Figure 6.19 Initial area decomposition that indicates the starting point in the top part of the waypoints and the trajectory. 6.2.1 Sequence reconfiguration simulations In the prototyping phase, the function was implemented in a similar way than the wind estimation and feature identification algorithms. The sequence reconfiguration consist in the redefinition of the order of the waypoints of a given surveillance mission to maximize the wind energy. This section describes the validation test performed to this functionality. Fort this test, the greater Seattle area was considered (see Figure 6.18). Since the area is considerably big, it was partitioned and divided into task that have to be performed by different UAV’s. For this case, the area was divided considering two mission agents and the resultant waypoint sequence of the first agent can be seen in Figure 6.19. Once the SITL simulation is performed, the summary of the resulting mission is observed in Table 6.8, in which no wind information was fed into the reconfiguration function prior 118 Chapter 6. Systems validation and verification, real experiments and results analysis Latitude (deg) 44.939 44.94 44.941 44.942 44.943 44.944 44.945 44.946 44.947 44.948 44.949 Longitude (deg) 14.49 14.495 14.5 14.505 14.51 Figure 6.27 A close view of the center area of Figure 6.26, the different curves depicting the effects on tension changes, which are very significant only in the case of sharp turns. Latitude (deg) 44.939 44.94 44.941 44.942 44.943 44.944 44.945 Longitude (deg) 14.484 14.486 14.488 14.49 14.492 14.494 14.496 59 61 42 60 68 23 41 67 22 66 40 21 39 Figure 6.28 A close view of the center area of Figure 6.26 when the desired curve considers the alignment of the wind to the slope unitary vector. This could be considered as the base of the dynamic soaring, since the Rayleigh cycles will require attitude changes. greatest concentration of points, it can be the effect on the interpolation splines can be observed more in detail while considering the effect of changing the slope to the direction of the next waypoint. This is depicted in Figure 6.28. If some of the waypoints are aligned to a theoretical wind direction, a close up effect on the resultant trajectory of the middle area can be seen in Figure 6.27. From the experiment performed in , a first solution on an easy implementation of continuous Rayleigh cycles among mission waypoints was extracted. If the resultant horizontal profiles have a sinusoidal vertical profile, then the resultant curve may produce a Rayleigh cycle. The next step was the consideration of the generation of the trajectory spline and the vector-field based tracking considering the effects of the interaction of the ocean and the atmosphere in the boundary layer. For this, a fully developed ocean was considered; the 6.2 Trajectory-related simulations 119 wind direction X Position (m) Y Position (m) Z Position (m) (a) wind direction desired course X Position (m) Y Position (m) Z Position (m) (b) Figure 6.29 Two test cases considered in the simulation of the full Rayleigh cycles. In (a), 2 full Rayleigh cycles are performed: starting from the dive-out maneuver at full throttle, the aircraft aims to dive-in into the next wave below the leeeddy formation zone to maximize the energy gain, the tension τi was adjusted manually at every cycle to match the wave frequency. In (b), 3 Rayleigh cycles are performed, the actual dive-in/dive-out maneuvers are winder with a general diminution of tension value τi compared to the first case. Note that in both cases the cycles are performed in a cross-wind direction. main objective was to assess the energy gains and losses that occur after performing 2 or 3 Rayleigh cycles. The two cases are illustrated in Figure 6.29. Once the simulations are performed, the next step is to analyze the airspeed results, which are presented in Figure 6.30. With the results presented in Figure 6.30 at this stage it is very difficult to conclude wether the changes on airspeed will lead to an overall energy gain, specially in the context of a full mission. However, the impact on changing the curve parameters is significant. In addition, it was concluded that the air/ocean interaction models have to be improved, since the observations performed by Richardson [ 23 , 24 ] suggest that the airspeed changes must be more dramatic, which will be translated into an actual energy gain. 120 Chapter 6. Systems validation and verification, real experiments and results analysis Step 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 Airspeed (m/s) 14 16 18 20 22 24 26 (a) Step 0 2000 4000 6000 8000 10000 12000 Airspeed (m/s) 12 14 16 18 20 22 24 26 28 30 (b) Figure 6.30 Airspeed of the UAS after the execution of the Rayleigh cycles depicted in Figure 6.29. (a) shows smaller airspeed changes, compared to (b), in which a lower tension value leads to greater changes in airspeed. 6.3 Validation Experiments The final step in the development of the FDES was an experimental campaign of several months that assessed the different functionalities on-board aiming to collect a big amount of data in different places and conditions for this analysis. This section will describe the details on such experimental campaign, which was performed in the southern region of Spain and in different seasons of the year, in which every aspect of the system was tested to evaluate its behavior and the proper integration with the system. 6.3.1 Experimental platforms There are two experimental platforms that were utilized in the various experiments of this test campaign. 6.3 Validation Experiments 121 Figure 6.31 The SkySurfer X8 is a 2000 mm wingspan fixed-wing gliding UAS; manufactured with EPO, it has a maximum takeoff weight of 1.7 kg. The fuselage length is 1.35 mm, it has a total surface area of approximately 46.3 dm 2 and a loading area of approximately 29.2 dm 2 , its center of gravity is located at 75 mm from the leading edge. The first one was The SkySurfer X8 (see Figure 6.31), which is a gliding UAS, mainly used for data collection for off board data processing but also for the online test of the wind field estimation and also the wind feature characterization algorithms. . The SkySurfer X8 was equipped with a PX4 controller, the actuators of the control surfaces are 6 servos of 17g, with a torque of 1.8kg/cm, and a speed capable of turning 60 º in 0.18s. The utilized battery was an 11.1v 2200mAh (20C) Lithium-Polymer battery, which powers the autopilot, the actuator and a brushless motor 3536-KV1200. Due to its limited space and payload the wind estimation function, when was performed online, the algorithms were calculated at the ground station, taking the instantaneous information from the telemetry. The second platform is the Proskytec V1 (see Figure 6.32), this platform is a V-Tail platform with better capabilities of both, weight and endurance, allowing testing at various levels of the different functionalities of the FDES. The Proskytec V1 was equipped with a PX4 controller, with two 9 g servos and two 17 g servos. It was equipped with a 22.2v 10000 mAh battery for maximum endurance. Also, the FDES module, the Odroid XU4, was also packed into the avionics compartment. This platform result ideal for the full implementation of the different functionalities and due to its relatively large payload, allows additional mission-related payload without any drawbacks related to the CG or to the takeoff weight. 122 Chapter 6. Systems validation and verification, real experiments and results analysis Figure 6.32 The Proskytec-V1 is a 2000 mm wingspan fixed-wing V-tail UAV; manufactured with EPO, it has a maximum takeoff weight of 3.5 kg, also it has a fuselage length of 1220 mm, a surface area of 70 dm 2 . Its center of gravity is located at 100mm from the leading edge. (a) (b) Figure 6.33 An aerial view of the Brenes facility (b), and the SkySurfer X8 (b) in a flight on this area. It can be observed that there are no bigger obstacles which are ideal for the analysis of the wind in a wide extension. 6.3.2 Experimental sites Two main sites were chosen as the testing facilities for the FDES modules, the first one is located in the Northern Metropolitan Area of Seville (37 ° 31’08.6"N 5 ° 51’31.5"W), close to the city of Brenes. This testing site has a wide area without any obstacles, such as mountains or threes, which could affect significantly the surface roughness, impacting the wind estimation. An aerial picture of this site can be observed in Figure 6.33. The second testing location, was the beach of El Terrón, located at the south of Huelva, 6.3 Validation Experiments 123 (a) (b) Figure 6.34 Images of the experiments performed in the region of Huelva, with the Proskytec V1 flying in the horizon while performing a soaring maneuver. Andalusia, Spain (37 ° 12’26.0"N 7 ° 10’24.5"W), which is a beach without population, which was ideal for the coastal tests. The winds at this region are predominantly coming from the west, and the average wind is of around 10 knots with typical gusts that can reach the 30 knots during the winter. Pictures of the tests conducted in this region can be observed in Figure 6.34. 6.3.3 Tuning experiments The PX4 autopilot, offers a control system based on a standard control with two main loops. The outer loop, also known as the guidance loop aims to determine the positional error with respect the waypoint in a mission (under the mission mode), and the inner loop which defines the control commands that are going to be given to the different actuators to move the surfaces. The architecture of such system is depicted in Figure 6.35. The PX4 manual suggests a general parameter configuration, based on the work presented by Meier et. al [ 83 ] and Park et. al [ 84 ]. The Park method uses a L1 control [ 77 ], and it has proved to be relatively efficient since the guidance logic uses the vehicle’s inertial speed for determining the acceleration command, which helps on the rejection of some disturbances, such as wind. After some experiments in the first testing site, it was determined that most of the suggested values in [ 75 ] are adequate for low winds. However, they have room for improvement, since the low winds ( µW<=10 knots) allow the controller to have a more aggressive response without unwanted oscillations. The baseline calibration was determined with the suggested values of the PX4 tuning guide [ 75 ] and a temperature of 25oC . The goal was to modify different parameters to achieve the optimal characteristics of the flight. These parameters are: the feed-forward roll gain Kφ f f , the proportional (roll-rate) gain Kφ P , 124 Chapter 6. Systems validation and verification, real experiments and results analysis Figure 6.35 The flight control strategy of the PX4 shows the two loops of control. This architecture was presented in [ 75 ], and shows the different stages for the different control modes that are available. the integral (roll trim offset) gain Kφ I , the feed-forward pitch gain Kθ f f , the proportional (pitch rate) gain Kθ P , the integral (pitch trim offset) gain Kθ I , and the acceptance radius Ωa . However, during the experiments, it was found that the reduction of the acceptance radius does not provide a significant impact on the performance, compared to the impact of changing the L1-period τL1. The resultant parameters were adjusted according to the increases suggested in [ 75 ], until the desired behavior was achieved. Table 6.10 show the baseline parameters, and the enhanced parameters for the given conditions. The selected mission for tuning was designed in conjunction with other fellows of the MarineUAS project. The primary goal was to solve a coverage problem; nevertheless, the type of mission suited perfectly the objectives of this research. The task consisted in visiting 55 waypoints, intending to coverage the green portion of Figure 6.33a. This mission is shown in Figure 6.36. After performing the different experiments, the tracking of the different waypoint showed significant differences with the selected values. The system was able to reject the wind disturbances, and the positional error was mainly related to the tuning parameters (see Figure 6.37). 6.3 Validation Experiments 125 Table 6.10 Tuned Parameters for the SkySurfer and Proskytec V1 in the region of Brenes. Parameter Baseline Value Enhanced Value Kφ FF 0,20 rad/s 0,36 rad/s Kφ P0,060 rad/s 0,075 rad/s Kφ I0,010 rad 0,031 rad Kθ FF 0,32 rad/s 0,44 rad/s Kθ P0,032 rad/s 0,46 rad/s Kθ P25 m 16 m Ωa2.0 m 0.1 m Figure 6.36 The selected mission that was utilized for the different tuning experiments. This mission was executed with the SkySurfer X8 UAS. The wind magnitude estimates, calculated with the direct computation method during the tuning experiments can be observed in Figure 6.38. It was observed that there is an initial jump on the wind speed magnitude. This is attributed to the convergence period of the estimation of the airspeed components at the autopilot. However, the average of the wind speed is relatively close to the average speed reported at the Seville Airport during the days of the flight, which were executed on November 6 to November 10 of 2017. During the 3 days, the average wind velocity was reported to be between 18 knots and 25 knots, with gusts that reached up to 35 knots, in accordance to the Sevilla Aeropuerto station. These experiments were performed at an altitude of 50 meters above ground, or 126 Chapter 6. Systems validation and verification, real experiments and results analysis Figure 6.37 The different trajectories and positional errors with the different tracking parameters in 4 flights. The black line represents the aircraft trajectory, the red area represents the error with the low parameters with a higher airspeed, and the green line is the error with higher airspeed with the optimized tuning parameters. approximately 80 meters. There are different jumps which can be attributed to the presence of discrete gusts among the flight. The FDES detected these as discrete gusts; furthermore, the hops may be attributed also to the changes in heading due to the sharp turns, that cause a sudden increase due to a poor definition of the angle limits on the PX4 native software. The heading of one of the experiments which is similar to the others can be observed in Figure 6.39. It is easily observed that some of the jumps on the heading calculation are directly cascaded into the wind estimation. Therefore, a correction needs to be performed to avoid this error and to determine the actual wind magnitude. This sudden jumps of the heading are attributed to the assumption of no sideslip in the calculation of the wind, and the direct use of the sin and cosine functions in the calculation, without considering the range limit changes. The next step was to produce a corrected wind speed based on the observations made at the previous stage. This correction was performed only by changing the input on the airspeed estimation of the autopilot by modifying the portion of the code that was illustrated in Code 2.1, which consider a zero sideslip into an approximation of the sideslip angle by comparing the ground speed/track vector to the airspeed/heading vector. The corrected wind speeds are observed in Figure 6.40 6.3 Validation Experiments 127 0 100 200 300 400 Time (sec) 0 5 10 15 20 Wind Speed (m/s) (a) Average wind speed of 11.7 m/s. 0 100 200 300 400 Time (sec) 0 10 20 30 Wind Speed (m/s) (b) Average wind speed of 12.3 m/s. 0 200 400 600 800 Time (sec) 0 5 10 15 20 25 Wind Speed (m/s) (c) Average wind speed of 13.6 m/s. 0 200 400 600 Time (sec) 0 5 10 15 Wind Speed (m/s) (d) Average wind speed of 10.7 m/s. Figure 6.38 This picture shows the estimation of the wind speed with the direct computation method for the four flights, with its trajectories were depicted in Figure 6.37. 0 200 400 Time (sec) -200 -100 0 100 200 heading (deg) Figure 6.39 The heading vs. time plot of the first run of the tuning experiments, the jumps of the heading are due to the crossing the 180ºand -180ªlimit. The corrected estimations were subject to a Weibull distribution analysis, Figure 6.40, from which it was determined the results that the estimates at the mission altitude are indeed Weibull-distributed data. And since the results show consistency with the weather