Intelligent Online Water Quality Monitoring through multiple Autonomous Surface Vehicles using the Bayesian Optimization Framework

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 Centralized Coordinator Bayesian Optimization-basedServer Autonomous Surface Vehicle1 Autonomous Surface Vehicle2 Autonomous Surface VehicleV Locations where to perform measurents Measurement Values, Location Information Multiple Surrogate Water Quality Parameter Gaussian Process-based Models Tesis Doctoral Doctorado Internacional en Ingeniería Electrónica Intelligent Online Water Quality Monitoring through multiple Autonomous Surface Vehicles using the Bayesian Optimization Framework Autor: Federico Daniel Peralta Samaniego Directores: Daniel Gutiérrez Reina Sergio Luis Toral Marín Ingeniería Electrónica Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2022 Tesis Doctoral Doctorado Internacional en Ingeniería Electrónica Intelligent Online Water Quality Monitoring through multiple Autonomous Surface Vehicles using the Bayesian Optimization Framework Autor: Federico Daniel Peralta Samaniego Directores: Daniel Gutiérrez Reina Profesor Titular Sergio Luis Toral Marín Profesor Catedrático Ingeniería Electrónica Escuela Técnica Superior de Ingeniería Universidad de Sevilla 2022 Tesis Doctoral: Intelligent Online Water Quality Monitoring through multiple Autonomous Surface Vehicles using the Bayesian Optimization Framework Autor: Federico Daniel Peralta Samaniego Directores: Daniel Gutiérrez Reina Sergio Luis Toral Marín 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 Fecha: A las personas que no dudan en hacer un esfuerzo para mejorar al mundo Agradecimientos Escribir un libro de tesis doctoral no es una acción individual, pero el esfuerzo de un grupo de gente. De igual manera, mejorar el mundo también lo es. Agradezco a mi familia, en especial a mi padre y madre, por enseñarme esto. Agradezco a mis profesores Daniel Gutierrez Reina, Sergio Toral Marín, por demostrarme que esto es cierto, además de agradecerlos por la enseñanza academica y profesional que fue impartida por ellos desde el primer momento de haberlos conocido. Agradezco a mis amigos y compañeros, por motivarme a seguir creyendo que mejorar el mundo es un esfuerzo colectivo, en donde todos tenemos nuestro espacio y nuestras oportunidades de seguir evolucionando como civilización, por y para el mundo. Federico Daniel Peralta Samaniego Sevilla, 2022 III Short Contents Resumen VII Abstract IX Short Contents XI Acronyms XVII Notation XIX 1 Introduction 1 1.1 Background and Motivation 1 1.2 Scope 6 1.3 Objectives 6 1.4 Thesis Contributions 7 1.5 Document Organization 8 2 Literature Review 9 2.1 Autonomous Vehicles 9 2.2 Global Path Planning 14 2.3 Modeling and monitoring applications 15 3 Methodology 21 3.1 Preliminaries 22 3.2 Bayesian Optimization Framework for Monitoring 29 3.3 Gaussian Processes as Surrogate Models 30 3.4 Acquisition Functions 37 3.5 Multi-Water Quality Parameter Estimation 42 3.6 Multiple Autonomous Surface Vehicles 45 3.7 Proposed Multi-ASV Multi WQP monitoring system 47 4 Implementation 51 XI XII Short Contents 4.1 Map Model 52 4.2 Water Quality Parameter Model Maps 52 4.3 Autonomous Surface Vehicles 54 4.4 Implementation using Python 57 5 Results 59 5.1 Simulation Setup 59 5.2 Kernel Selection 60 5.3 Acquisition Function Selection 63 5.4 Multi-ASV 69 5.5 Multi WQP monitoring using Multi-ASVs 69 5.6 Comparison with other methods 74 5.7 Summary of the Results 78 6 Conclusion and Future Work 81 6.1 Conclusion 81 6.2 Future Work 83 6.3 Publication List 84 Bibliography 85 Contents Resumen VII Abstract IX Short Contents XI Acronyms XVII Notation XIX 1 Introduction 1 1.1 Background and Motivation 1 1.2 Scope 6 1.3 Objectives 6 1.4 Thesis Contributions 7 1.5 Document Organization 8 2 Literature Review 9 2.1 Autonomous Vehicles 9 2.1.1 Autonomous Systems 10 2.1.2 Autonomy 12 2.2 Global Path Planning 14 2.3 Modeling and monitoring applications 15 3 Methodology 21 3.1 Preliminaries 22 3.1.1 Map model and measurement locations. 22 3.1.2 Water Quality Parameter Model Maps. Measurement values 24 3.1.3 Autonomous Surface Vehicles: Agents and coordination 25 Guidance, Navigation & Control (GNC) 25 ASV constraints 25 Centralized Coordination 26 Water Quality Parameter Sensor System 26 XIII XIV Contents 3.1.4 Objective Functions 27 Performance Metrics 28 Multi Objective Optimization 29 3.2 Bayesian Optimization Framework for Monitoring 29 3.2.1 Surrogate Model 30 3.2.2 Utility Function 30 3.3 Gaussian Processes as Surrogate Models 30 3.3.1 Covariance Functions 31 Example of covariance functions 32 Operations Between Kernels 34 3.3.2 Solution of a Gaussian Process 35 3.3.3 Models and Hyperparameters Selection 36 3.4 Acquisition Functions 37 3.4.1 Proposed Acquisition Functions 40 Adaptation functions 40 Definition of maximum distance for adapted measuring 41 3.5 Multi-Water Quality Parameter Estimation 42 3.5.1 Multi-objective problem 43 3.5.2 Multi-function adaptation 45 3.6 Multiple Autonomous Surface Vehicles 45 3.6.1 Initial Vehicle Positioning 46 3.6.2 Region Partitioning using Voronoi Diagrams 46 3.7 Proposed Multi-ASV Multi WQP monitoring system 47 4 Implementation 51 4.1 Map Model 52 4.2 Water Quality Parameter Model Maps 52 4.2.1 Benchmark Functions 53 4.3 Autonomous Surface Vehicles 54 4.3.1 Guidance, Navigation & Control (GNC) 54 4.3.2 ASV Constraints 56 4.3.3 Water Quality Paremeters Sensor System 56 4.3.4 Centralized Coordinator Server 57 4.4 Implementation using Python 57 5 Results 59 5.1 Simulation Setup 59 Simulated Water Quality Model using Benchmark Functions 59 Theoretical maximum energy available 59 5.2 Kernel Selection 60 5.3 Acquisition Function Selection 63 5.4 Multi-ASV 69 5.5 Multi WQP monitoring using Multi-ASVs 69 Contents XV 5.6 Comparison with other methods 74 5.7 Summary of the Results 78 6 Conclusion and Future Work 81 6.1 Conclusion 81 6.2 Future Work 83 6.3 Publication List 84 Bibliography 85 Notation Acronyms AAV Autonomous Aerial Vehicle AF Acquisition Function AGV Autonomous Ground Vehicle ASV Autonomous Surface Vehicle AUV Autonomous Underwater Vehicle AV Autonomous Vehicle BO Bayesian optimization CC Centralized Coordinator DRL Deep Reinforcement Learning EI Expected Improvement GA Genetic Algorithm GNC Guidance, Navigation & Control GP Gaussian Process GPP Global Path Planning GQI Good Quality Indicator IPP Informative Path Planning LPP Local Path Planning MFE Multi-Function Estimation MiP Mission Planning MOO Multi-Objective Optimization MRS Multi-Robot System NSGA-II Non-dominated Sorting Genetic Algorithm II PESMOC Predictive Entropy Search for Multi-Objective BO with Constraints PI Probability of Improvement POMDP Partially Observable Markov Decision Process PP Path Planning XVII XVIII Chapter 0. Notation PSD Positive Semidefinite PSO Particle Swarm Optimization RBF Radial Basis Function RQ Rational Quadratic RRT* Rapidly-Exploring Random Tree star SE Squared Exponential SEI Scaled Expected Improvement SF Shekel function TG Trajectory Generation TSP Traveling Salesman Problem WQP Water Quality Parameter XIX Notation RnReal Number Hyper-space MMatrix Model Map σ(x)Standard Deviation value of a GP Regression p(y|x)probability of ygiven x α(·)Acquisition Function ` Length-scale hyperaparameter for SE, RQ and Matérn Kernels εUnknown noise value exp·Exponential function Γ(ν)Gamma Function λProportion of length-scale value log Logarithmic function EExpected value V[x]Normalizing function for SEI AF µS/cm microSiemens per centimeter µ(x)Mean value of a GP Regression Mi,jElement of the map model N(·,·)Normal distribution function XFeasible Search Space φ(x)Probabilistic Density Function Φ(Z)Cumulative Distribution Function σ2 nExpected measurement noise variance ~vVelocity of an ASV ξExploration bias hyperparameter posterior subscript for posterior component of a GP prior ,psubscript for prior component of a GP ||·|| Norm function LLikelihood of a distribution ALocation of peaks in a SF cInverse powers of a peak in a SF dside length of Mi,j DDistance traveled by an ASV d(·)Euclidean distance between different x eAverage current draw for ASV movement f(x)Unknown target Function F(x)Objective Function ˆ f(x)Approximation of target function ¯ fs(x)Mean value of the real data xBOptimized (Maximum/Minimum) measurement location g1(x)Feasibility check constraint function IIdentity Matrix KKernel Matrix k(x,x0)kernel evaluation 6Chapter 1. Introduction work because it defines the movements of the ASV. Thus, classical acquisition functions used in BO cannot be directly applied to the monitoring problem due to the mobility restrictions of an ASV. In this work, we proposed several adaptations of the classical acquisition functions that are more suitable for monitoring tasks. To undergo the process of eutrophication, water bodies must suffer from several atypical changes. These values can be obtained beforehand (or during eutrophication) in order to prevent (or recover from) water pollution and prevent algae blooms. In this sense, multiple WQP values need to be obtained from the current state of water bodies. Therefore, an important factor is the simultaneous measurement of several WQPs, so an efficient multi-function estimation system is proposed, as well as a Multi-Objective Optimization (MOO) strategy solved with GA. It is used to combine different utility functions because it is an optimal decision-making strategy whenever there are several objectives, as in this case study. Furthermore, in this work, we make use of different tools and techniques to make monitoring with an ASV possible, such as constraint handling and adaptations of AFs. Note that data acquisition should be cooperative, i.e., several ASVs are available and need to jointly obtain measurements. Therefore, an intelligent multi-ASV system is proposed, where the ASVs are in charge of efficiently performing measurements using a Voronoi Partition system to effectively create online coverage regions. The Voronoi Partition or Voronoi Diagram (VD) is a mathematical expression that defines regions according to a set of generators [ 21 ]. If the positions of the ASVs serve as generators, each region has an ASV assigned to it, and the BO technique can be utilized seamlessly for each ASV to choose new measurement locations within each of the defined regions. 1.2 Scope Monitoring WQPs can be interpreted as a recurrent task aimed at obtaining useful information about the target environment. In this thesis, a single mission is presented, as it will produce the current state of various WQPs of a water body using several ASVs for performing measurements. All ASVs are supposedly equipped with the same systems, which implies that these agents are homogeneous and will present the same behavior and constraints. For the selection of the different components, several experiments will be performed so that the system proposed in this thesis and its results are replicable and applicable to any water body. 1.3 Objectives The previous sections served to summarize the background, general task, vehicle usage, and scope of this thesis. The General and Specific Objectives of this thesis are described below, with keywords highlighted in bold. The General Objective is as follows. Design and implementation of a multi-autonomous surface vehicle monitoring system aimed at obtaining models of multiple-water quality parameters using a Bayesian optimization framework with a Voronoi partitioning system. 1.4 Thesis Contributions 7 The specific objectives are designed to separate the general objective into specific components that the system needs to ensure efficacy, efficiency, and robustness. They are are described as follows. • Propose a Bayesian optimization-based Monitoring Mission for a multi-vehicle system in unknown aquatic environments. • Perform intelligent autonomous surface vehicle deployments in order to obtain a set of real or simulated approximated water quality parameter maps of a water body using Gaussian Processes as Surrogate Models, using a fleet of surface vehicles. •Evaluate and compare the proposed method with common exploration/monitoring methods. 1.4 Thesis Contributions The main contributions of this thesis can be summarized as follows. • Application of Bayesian optimization for sequential decision making in a Water Quality Monitoring scenario, where the general BO equation is used to obtain measurement locations that will be useful for generating WQP models according to Gaussian Process stochastic models based on kernels or covariance functions. • Hyperparameterization and comparison of classical kernel functions as main components of stochastic surrogate models for water quality distribution maps, validated through the usage of expected WQP behaviors in simulations using benchmark functions. • Proposition of three adaptations of acquisition functions based on classical techniques for monitoring missions using an ASV with energy constraints, so that data can be acquired taking into account the limitations of mobility of the ASVs in charge of monitoring. • A centralized Voronoi-based coordination system for multiple ASVs based on active region partitioning and data sharing, to create robust operations, seamless cooperation, and achieve better results in less time. • A generalized multi-parameter measuring system for environmental monitoring based on Bayesian optimization with fusions of functions, considering centralized communication of the different ASVs, asynchronously. • A multi-model acquisition system through multiple ASVs for online environmental monitoring. The system considers not only the existence of muti-WQPs but also multi-ASVs. • Performance evaluation of the system considering the coefficient of determination of the obtained models with benchmark functions. Comparison of the proposed system with related monitoring works. 8Chapter 1. Introduction 1.5 Document Organization The rest of this book is organized according to a series of chapters that contain the following topics: Chapter 2 presents the related work and starts by defining the basic concepts of real components of the system according to the current state of the art, including the definitions of autonomous surface vehicles. This chapter describes the usage of other mathematical components of the system, such as Bayesian optimization and multi-objective optimization. Finally, current efforts for Monitoring Applications are also described in the chapter, addressing not only general water quality monitoring, but also monitoring research and implementation for highly eutrophized water like Ypacarai Lake. Next, Chapter 3 describes the main architecture and each part and component of the system. This includes considerations, assumptions, and constraints according to the scope defined previously. Bayesian optimization is properly described, as well as its main components: Gaussian Processes and Acquisition Functions, their mathematical definitions, limitations, and solving methods/strategies. It includes examples of covariance functions used in the work. Voronoi regions are described in this chapter, as well as Multi- Function Estimation and Multi-Objective Optimization strategies. These components are fully described in this chapter, but their implementation is described in Chapter 4, where the usage of numerous pre-built systems is explained. The embbeded systems that implement the system are also described in this chapter. The results are shown in Chapter 5. It includes complete results in detail by component, as well as summarized results. Several tables can be found that show the numerical performance values obtained, as well as a discussion of the results section. Comparison with other methods has also been done, the implementation and results are also found in this chapter, accompanied by a discussion of the obtained results. Finally, Chapter 6 concludes the work showing the learned insights, the overall results of the system, and future work. 2 Literature Review Galaxies make stars, stars make worlds. Neil deGrasse Tyson, 2020 2.1 Autonomous Vehicles Ev ery year, more robots perform their tasks and objectives with fewer margins of error; this is the result of countless research and advances in each of the components that generate what is known as an autonomous system. These components can be grouped into three main systems according to Liu et al. [ 22 ]: Navigation, Guidance and Control systems, known as GNC systems. These components answer some critical questions [ 23 ]: “Where am I?” , “Where do I want to go and how?” , and “What can I do to get there?” . Robots that can answer these questions can act on their own in order to accomplish their goals, and therefore be perceived to have a higher level of autonomy. Strategically, AVs research can be done aiming to answer these questions and develop systems that provide a solution to theses problems, specifically the ones that involve navigation, guidance and control. Searching for a solution to the latter issue, Control Systems have been in development since much before the birth of robotics, these control systems are used to determine appropriate control forces to be generated by the actuators, according to Liu et al. [ 22 ], and are already very well designed because of its necessity as a way to control any machine, including autonomous ones. Control Systems are in charge of the movement, which is a fairly well studied subject. A common situation is that these machines are in the need of constantly identifying their states and the states of their surroundings, this identification is accomplished by a Navigation System, which uses sensors and mathematical models to perform an observation of the situation and/or calculate the position of itself in respect of the environment (e.g. answers the question "Where am I?".) Research on movement and observation has been extensively developed in this field, but planning is a relatively new area and is the missing piece in the Autonomy 9 10 Chapter 2. Literature Review puzzle [ 22 ]. Planning is what Guidance Systems are responsible for, and even though in the past, the useful signals produced by Navigation Systems and the ones that Control Systems receive did not have a computer as an intermediary, nowadays these systems work and cooperate seamlessly due to the implementation and use of Autonomous Systems, which are processing and producing the necessary signals. Clearly, current research on AVs is mainly focused on developing optimal deployment, intelligent motion, or smart guidance system, these terms correspond to planning [24, 25, 16]. 2.1.1 Autonomous Systems The deployment of vehicles are always mission-dependent; that is, when the deployment performance is defined, it always refers to the level of compliance of the mission with respect to the actual movements of the vehicles. This premise limits the vision of the average researcher to the development of Path Planning (PP) subsystems (or components). On the other hand, most of high performance Autonomous Systems also manage to include an appropriate Mission Planning (MiP) component that modifies the way of searching next goals or locations. Another way to improve the global effectiveness of the autonomous system is to focus on Trajectory Generation (TG) after the PP procedure. TG ensures that the control system positions the AV in known locations. Based on the work by Liu et al. [ 22 ], Autonomous Systems can be defined with a set of the 3 components introduced in this paragraph (Figure 2.1.) It is also worth noting that the system has an input signal “AV and Environment Information” and an output signal “Desired Velocity”, both of which were previously noted as output and input from the other systems that comprise a GNC or Autonomous System. When it comes to the problem of developing GNC systems, one of the first decisions involves choosing the right method to solve global optimization problems. Once the domain of the environment is defined, the system must find the global minimum of a function f(x) , which can be a defined function or a function not analytically representable. In either way, a solution must be found and there are two general types of approaches: stochastic and deterministic [ 20 ]. Whenever research papers present a new way of solving a particular problem or a review of methods and techniques found in the state of the art, it is usual for the algorithms to be compared with algorithms of the same type (stochastic or deterministic), and this can be considered a good practice because it is assumed that the other approach may underperform with respect to the objective of the research. It is known that stochastic methods will identify the global minimum of a function in an infinite time with probability of one, and that deterministic methods ensure that with exact computations and arbitrary long times, an approximation of the global minimum will be found after a finite time. All in all, Guidance Systems are designed and implemented, with a focus on the mission of the project under development. It is necessary to define and exemplify the components of the GSs in order to deeply comprehend the Guidance System design paradigm and, therefore, propose and produce better guidance systems. • Mission Planning (MiP): referring to global high-level decisions that need to be made, Liu et al. [ 22 ] define that a MiP component should obtain all possible data (tasks, environment, operators, etc.) and define missions accordingly. These missions can 2.1 Autonomous Vehicles 11 Figure 2.1 Autonomous System components: Mission Planning, Path Planning and Trajectory Generation. be to obtain a target point or proceed to perform a task. Zhen et al. [ 26 ] describe MiP as a Waypoint (or target) Generation module for a multi-UAV system, Sampedro et al. [ 27 ] differ in these definitions, defining a Global Mission Planner (which acts like a normal MiP by providing waypoints) and an Agent Mission Planner that is in charge of executing the tasks made by the previous planner. In the survey by Mac et al. [ 28 ], mission planning can arguably be defined as the “cognition” step in the autonomous process they define. MiP component is among the more abstract components of AVs not only because they need to provide context to the lower levels of Autonomous Systems but also they need to transform the abstract ideas of Mission and Vision Statements into feasible actions. • Path Planning: The process of obtaining a valid route between two points is called Path Planning and is the crucial component of the GS due to the fact that to safely guide one vehicle to a point, a route must be planned. Zhou et al. [ 23 ] define this stage as the route planning procedure and also as a macro-modality constraint problem, and the vehicle is treated as a particle in the environment. Path Planning takes in all the data from NS and the MiP component, and obtains a series of waypoints to visit in order to reach to a desired location. Shin and Chae [ 29 ] elaborate a performance review of some of the best known path planning algorithms, such as A* or GA, using 12 Chapter 2. Literature Review a small grid-based experimental configuration and a simulated discrete particle as vehicle. Other works elaborate similar performance reviews but take into account models of the vehicle dynamics to test the algorithms, for example, Chaari et al. [ 30 ] using a simulated Autonomous Ground Vehicle (AGV) while Zammit and Van Kampen [ 31 ] used a simulated AAV as vehicle for comparison. These works define various types of path planning methods, they can be Global Path Planning (GPP) or Local Path Planning (LPP), the first approach takes into account only information found in the past (be it from the environment, or even from the vehicle itself), it can be considered as a first step towards a high-performance planning or a previous step for LPP, which is the approach that considers continuously updated information and obtains routes that can change through time if required. Though Path Planning is sufficient to produce valid routes for vehicles, in some occasions additional computation is needed in order to provide safe routes. • Trajectory Generation: Whenever a more complex scenario is presented, this component is necessary to ensure that the vehicle follows the planned or projected path. More complex scenario means more complex dynamics and/or more rules to be followed to reach a destination. Panden et al. [ 32 ] denote that the TG component is suitable whenever the environment or constraints are dynamic, this is true in most real-life scenarios since environment hazards such as winds, waves, etc. are present. In spite of the vast amount of definition for these systems, all researchers in this field define these systems in the same way as stated in this thesis, differing only in word selection. For the sake of clarity, in this thesis, the terms from different works will be adapted to the previous definition whenever necessary. 2.1.2 Autonomy Figure 2.2 is designed around the definition of autonomy, taking into account the components of the Guidance Systems. The less the AV decides its actions, the lower the autonomy level, therefore, the less it comprises all possible components of the Guidance Systems. However, the process shown in Figure 2.1 can be misleading because it shows that each component depends on previous ones, and while this is true, if a work attempts to work on more general components, it is assumed that more specific components are also developed, but this is not true. Some authors focus their work on a specific component like Path Planning and not Trajectory Generation, while others develop more components to ensure a higher level of autonomy. To emphasize the definition of autonomy, it is necessary to differentiate this term from automation, the first refers to the capacity to deliberate a course of action by itself, while the latter most often implies the process of performing physical or logical actions without any help. Zhou et al. [ 23 ] proposed eight levels of automation, in which each designed Guidance System can be assigned to each of the levels. However, it must be noted that these levels include cooperation between vehicles (i.e.: sharing information, etc.) from level 6+ and it is safe to say that a vehicle can be fully autonomous while not being able to cooperate with others. Abraham et al. [ 33 ] define seven levels of automation for AVs and the levels differ in definition but correlate in the desired output with the previous cited work. Wickens [ 34 ] 2.1 Autonomous Vehicles 13 Figure 2.2 Automation capabilities. made a third automation definition that is based on stages rather than levels and defines 4 global situations: Attention filtering, Diagnosis, Action Choice, and Action Execution. These definitions were consolidated in Table 2.1 to have the least possible number of levels of Guidance Systems but sufficient to differentiate the perceived quality of automation. Table 2.1 Automation Level. Automation General Level (Zhou Level (Abraham Stage Level Definition et al. [23]) et al. [33]) (Wickens [34]) 1 Human Guidance 1-2 1 1 2 LPP 3 2-4 2 3 GPP & LPP 4-5 5 - 4 PP & TG 6-7 6 3 5 MiP, PP & TG 8 7 4 Guidance Systems design must exploit the vehicle’s capabilities, properties and applications, i.e.: for maritime rescues purposes, the vehicle must prioritize time and power management, for unknown terrain mapping, the vehicle must cover the entire environment in detail, and for precision agriculture, the vehicle must carefully select the paths to drive. In Table 2.2, several examples of work in this field are shown. The approach that relates to this thesis is the work of Morere et al. [ 19 ], where an AAV deliberates the next goal position according to the level of uncertainty of the environment and aims to decrease the uncertainty in high-gradient zones. 14 Chapter 2. Literature Review Table 2.2 Automation level of some vehicles. Guidance System Perceived Automation Level Main Focus Vehicle Type Morere et al. [19] 4 MiP, PP, TG AAV Tokekar et al. [35] 4 GPP, LPP AGV & AAV Liu et al. [36] 3 LPP ASV 2.2 Global Path Planning Some applications require redundancy, dangerous activities, coverage, and emergence. These are defined as applications that can be achieved by a Multi-Robot System (MRS), which is more robust, flexible, and scalable than a single robot system (isolated vehicles) according to Bayindir [ 37 ]. MRS have a special place in the field of mobile robotics because they deal with the organization of many, organization that should be comprehensible from a general point of view. The last statement by no means implies a general supervisor or coordinator, MRS can be centralized or decentralized [ 37 ]. They are capable of providing a solution whether a central unit is providing the tasks or not. Many models of social robots and AVs are based on natural behaviors or the human interpretations of them, these models are focused on defining rules and instructions for a single vehicle, they can intentionally cooperate or compete in a social environment but in the first approach of development, then on a second stage of development, the rules are expanded and include interactions between vehicles in order to achieve the social requirements that are in the scope of the research. Addressing, for example, the coverage problem, it is common to have a vehicle that is capable to cover a pre-defined area and later develop the communication and cooperation system in order to divide the areas to cover when a MRS is present. Table 2.3 shows a number of examples of MRS, grouped by application. A coverage application aims to monitor a geographical area, swarming is the action of collaborating between simple robots to accomplish a goal, transport robots are vehicles that coherently act and communicate to move one large object. From the examples shown in Table 2.3, the only application that cannot be accomplished by a single vehicle is swarming, there are numerous works that aim to patrol, transport or even cover an area, but for evident reasons, they take more time and more prone to fail the task at hand. In any case, it is important to consider the research that fits the scope of this thesis, which includes all coverage application researches and works. The thesis by Kapoutsis [ 38 ] provides a multi-ASV approach to cover an unknown environment. It focuses on dividing the environment using a cyclic algorithm taking into account the starting position of the vehicles. While enclosed environments can be optimally covered in a finite time, open environments are difficult to cover because in some of the worst-case scenarios the supposed 2.3 Modeling and monitoring applications 15 Table 2.3 Examples of coverage applications. Application Vehicle Type Main Alg. References Coverage AGV Flocking [25] ASV DARP [38] Swarming AAV Forward, Coherence, Avoidance [39] Transport AGV Formation, Local and Global PP [40] previous information can lead to dead-ends or local minimum. This problem can be avoided using the Kapoutsis method proposed [ 38 ] or by continually updating the information about the environment, more specifically, updating the posterior model of the environment with the prior model and a new evidence taken from the environment. Although this method can be more computationally intense than the first one, it provides better coverage whenever the information of the environment is more important than the coverage of the environment itself. This thesis describes a method that aims to obtain approximate models of unknown environments in a shorter total time than current methods, admitting uncertainty. This will be accomplished by a group of Autonomous Surface Vehicles that distributively work to obtain information about the environment and mathematical models that take into account uncertainty and probability. The automation level will reach a level of 5 (Table 2.1), considering that this thesis addressed the Mission Planning component, and uses multiple systems to perform Path Planning and Trajectory Generation. This is one of the main contributions of this thesis, as initially described in the previous chapter. 2.3 Modeling and monitoring applications Regarding GPP for monitoring, there exist conventional monitoring methods that have been implemented due to their simplicity, like the lawnmower method. It consists of sweeping the space using a predefined distance between measurement locations (in the case of discrete measuring) and a preestablished route morphology. In general, using these techniques yield a slow but detailed monitoring. Therefore, these are the preferred methods when there are neither time-related restrictions nor expensive exploration costs. L-Cover, T-Cover, and Z-cover [ 16 ] are some common methods that fall into the lawnmower category. These methods are extremely useful in situations where the work or search space is very narrow, such as with rivers and streams; however, they have less impact observational wise as the water body widens due to energy availability concerns. In addition, these methods are not designed to do online learning or monitoring. Various coverage, patrolling, and 22 Chapter 3. Methodology Start End Obtain the current surrogate state of the environment Stopping Criteria met? Save new data, update current belief of the environment Evaluate the real state of the environment in the mentioned point Generate a new evaluation point according to internal/external states No Yes Figure 3.1 Abstract general flowchart of the steps that needs to be done in the monitoring mission. 3.1 Preliminaries The monitoring system designed in this thesis has the objective of efficiently obtaining water quality models of a given water body. A water body is defined as any significant accumulation of water on the surface of the Earth and can be a reference to seas, salt lagoons, freshwater lakes, rivers, or small lagoons. In this context, each of many water quality parameters is mathematically defined as an unknown black box target function f(x) , where the argument x corresponds to a location in the search space (water body). Evidently, each target function is unknown to the observer, and additionally, their true behaviors could never be obtained. Consequently, these functions are replaced by surrogate models designed to approximate the real behaviors. The values of the target functions are available to the system, but the system performs noisy reads such that: ˆ f(x) = f(x)+ε(3.1) with ε as the noise bounded to the measurement. With that, the expected output of the system will be based on these noisy measurements rather than on the true values. Additionally, for computation purposes, the real model will be replaced with a surrogate model that is set to behave like the real one. A comprehensive definition of the surrogate model is available in later sections, but, loosely speaking, the surrogate model can be thought of as a model of ˆ f(x)∀x∈X. 3.1.1 Map model and measurement locations. X is a subspace of the Rn space. It is expected that x is a geodetic coordinate, that is, a latitude, longitude, and altitude data, because the agents in charge of performing the measurements are on the surface of a water body. However, since the agents can only travel on the surface, the space of actions reduces to Latitude and Longitude. Moreover, 3.1 Preliminaries 23 since the dimensions of Earth will always be much larger than the dimensions of the water body to be monitored, the space X reduces to a 2D projection of the water body. Then, the modeled map can be described in more detail, since only the rectangle that encloses the water body is needed. This rectangle can be thought of as a matrix M , which can be called a map model. Therefore, the monitoring region is considered as a matrix M with dimensions of m×n , i.e, a monochrome image with pixel dimensions m×n , i.e., each pixel corresponds to a square of side d that has a black or white color that represents the occupancy state of the square. This square represents a compromise solution in which the environment is sufficiently detailed but large enough to allow large-scale modeling. Each element of the map model Mi,jhas a value defined by an occupancy state (1 if the element corresponds to a occupied square space d×din the space, 0 if not). Figure 3.2 Example of a map model. Originally, the real image is drawn from a Satellite Resource. Then, it is translated to a binary image according to the navigability inside a square of a defined side d. Figure 3.2 shows an example of the map model of a water body used in this work; in the image, the color white represents obstacles (including prohibited zones, natural obstacles, land, etc.), while the color black represents a square of free obstacles. Definitely, the environment studied is discrete due to the usage of the map model M , so that a vehicle can only be located in a position related to an element of the map. Finally, the set of navigable elements is defined as follows. X⊆M⊆R2(3.2) A general way to represent the map is to also include every non-navigable location in the X with the addition of a constraint function that yields the occupancy state. The constraint function will satisfy the following inequality for all feasible locations: g1(x)≤0 . Additionally, any location outside of the map model M will yield undefined values, i.e., the edges of the map model define the boundaries of the region in which an agent can be 24 Chapter 3. Methodology found. However, in what follows, unless defined otherwise, a generic location x∈X will be considered to satisfy the constraint. An important consideration is that the studied water environments are always non-convex sets. The convexity of a set can be described by graphical observations. If, in the space of actions (map model), there exists any pair of feasible locations that cannot be connected (through a segment line) without selecting unfeasible locations, the set has a non-convex property. Regarding the fact that any water body will generally have a non-standard shape, it is better to consider that the water body will always be a non-convex set, including the example shown in Figure 3.2. Consequently, local path planners are needed in this scenario to ensure that each available ASV can travel safely from one location to another. 3.1.2 Water Quality Parameter Model Maps. Measurement values Observing Eq. (3.1) , x , as well as f(x) , are now mathematically defined, but f(x) and the measurements ˆ f(x) must also be defined in terms of the monitoring application, as x was in the preceding paragraphs. Environmental water quality is evaluated according to a set of physicochemical and biological parameters of the water. In the case of water bodies, the parameters must meet a defined criterion so that they can harbor life and maintain a balanced healthy ecosystem. These parameters can be divided into three subclasses, physical, chemical, and biological indicators. Current available water quality sensors include both physical and chemical sensors, but not biological ones. Therefore, the available Water Quality Parameters or Indicators relevant to this thesis are included in Table 3.1. Note that the values for Good Quality Indicators are not defined for drinkable water but for freshwater [80] capable of sustaining life. Table 3.1 Water Quality Parameters (ordered by their importance according to [ 80 ]) available for measuring indicating their abbreviation, whether they are a chemical (C.) or physical (P.) indicator, unit of measurement, typical range, and typical Good Quality Indicator (GQI) value for freshwaters. Parameter Abvr. Type Unit Typical Range GQI Dissolved Oxygen DO C. % 0∼100 75 −125 Potential of Hydrogen pH C. H+ 0∼14 6 −8.5 Nitrate NO3 C. mg/L 0 - 15 ≤5 Temperature T P. ºC0∼100 — Turbidity Tu. P. NTU 0∼4,000 ≤20 Conductivity Cond. P. µS/cm 0∼60,000 200 −800 Oxidation-Reduction ORP C. mV 0∼2,000 300 −500 Potential A Water quality parameter model is defined using the function f(x) , a value or evaluation of the model at the location x is always available for measurement. The model represents the current behavior of a particular WQP and is considered continuous and multimodal. These characteristics are the expected characteristics of water quality parameter maps. 3.1 Preliminaries 25 WQPs are not necessarily time-invariant, implying that their stochastic behavior is not stationary since the environment is continuously changing. In that sense, a study on WQP modeling must also be done to obtain time frames in which the measurements are valid and can yield good results. 3.1.3 Autonomous Surface Vehicles: Agents and coordination The agents in charge of performing the measurements are ASV. The v th ASV, considering v the id number of an ASV in a fleet of V ASVs, i.e., v= [1,2,3, ... ,V] , must travel on the surface of a water body from location to location, loiter around the measurement point, i.e. stopping and managing to stand still, perform the measurement (and save/send/store), and obtain the next measurement location according to the system. All these actions must be performed by their on-board system. Therefore, it is necessary to define the specific on-board systems of the ASVs. Guidance, Navigation & Control (GNC) As discussed in the previous chapters, all ASVs are admittedly capable of traveling through the surface of any water body and, in this case, can also perform Water Quality measurements. The GNC system on board is capable of determining the location of an ASV can be determined by assuming a certain error. The Navigation subsystem provides a location value of the vehicle within the real world; however, since the positioning system is not perfect, the positioning error leads to guidance inaccuracies. These inaccuracies in the guidance include considering that an ASV is not at the estimated location. To efficiently overcome this issue, the map model to be obtained can account for this and define that the square side length d , previously defined, is greater than the expected error output derr of the location estimation. Then, the assumption is that the WQPs measured at the real location x are the same for all nearby locations that are d units away, in terms of latitude and/or longitude, including the assigned measurement location x∈M. The Guidance subsystem will plan paths according to the current location and the goal location. Taking into account the previous definition of this subsystem and its components, in this thesis, each ASV is considered to have a dedicated Local Path Planning component, which uses an already defined goal to obtain a collision-free path that can effectively and efficiently move the ASV from one location to another. If there are no obstacles in the segment, a direct route is planned; otherwise, the guidance system plans a path using RRT*, as it can provide quick and good routes. The control subsystem is the last subsystem of the GNC, and it actively modifies the motor and other actuators states to reach the desired locations. However, since the navigation subsystem has an error, the goal location may be different from the final location reached. This error also accounts for water and air currents, stopping mechanisms, etc. Consequently, the Mission Planning and Global Path planning components are the only non-essential components that are left to be defined, but will be in the next chapter. ASV constraints Regarding energy and autonomy, considering that the WQP sensor system is powered by its own battery and that the onboard electronics consume relatively very little power, battery usage, and time consumption are considered constraints, defining the maximum distance 26 Chapter 3. Methodology that each ASV can travel. The battery of the WQP sensor system can run for long sessions, noticeably longer than the expected mission time, while the onboard electronics draw current from the main battery, but at a significantly lower rate than the motors. Moreover, the mission planning and global path planning components do not run on each ASV, as they run on a server-side system. The usage of batteries per unit of movement cannot be exacted, but an approximation can be made considering the movement factors, the available battery, the motors utilized, their acceleration capabilities, and the desired velocity of the ASV. Considering that on average an ASV draws e current per unit of distance driven, the average velocity of the ASV is ||~v||, and the total distance driven is D, the energy used is: Tenergy =e ||~v|| ×D(3.3) Finally, with a defined maximum energy use Tenergymax , usually in terms of Amperes hours [Ah], the maximum distance Dmax can be obtained using the latter equation. This value could be used to express when the mission has ended. Normally, this maximum energy usage will be less than the energy available to ensure that an ASV can return to the shore of the water body after it ends its mission. Centralized Coordination The centralized coordinator (CC) is responsible for running the system. It stores the data obtained by each vehicle in a database and is responsible for providing an optimal measurement location for each vehicle. It is centralized because it serves as a global decision-making component and coordinates the movement of multiple available ASVs. The paradigm is centralized coordination and decentralized execution, which yields effective coordinate solutions (surrogate models that include information obtained by all vehicles) and efficient asynchronous measurement of water quality parameters through multiple vehicles. Water Quality Parameter Sensor System Each ASV has a system dedicated to measuring multiple water quality parameters (pH, DO, Tu., etc.). When an ASV reaches a desired measurement location, it stops and takes measurements of all available WQPs. There are S different sensors, each one is accessed through its id number s , considering s= [1,2,3,...,S] . To take measurements, each ASV must stop and instruct to perform a still measurement, preventing the socalled continuous monitoring, i.e., the vehicle cannot constantly obtain new information. Consequently, measurements are performed sequentially and, at the end of the mission, a discrete number M of measurements that have been performed is available. In that sense, pm, with m= [1,2,3,...,M], is a location x, where a measurement has been performed. Since multiple WQPs are measured simultaneously, it is necessary to describe multiple target functions fs(x) , and, consequently, multiple measurements ˆ fs(pm) at location pm . For simplicity, each measured WQP is referred to as fs(x) instead of their real name, e.g., the DO (Table 3.1) function for a water body is f1(x) , etc. Then, the mission of the proposed system pivots to selecting measurement locations pm so that each surrogate model ˆ fs(pm) is as accurate as possible to the real behavior, which in turn can be described in terms of a Multi-Objective Optimization. 3.1 Preliminaries 27 The specific required on-board system is described in Figure 3.3. It comprehends the 4 systems discussed in this section. The input values are directed to specific subsystems so that autonomy can be efficiently achieved. The output values are directed towards changing the location of the ASV or to inform the coordinator about measurements performed. In Figure 3.4, an overview of the behavior of the system can be observed. It is a simplified version that shows only what matters to the mission and objective of this thesis. In this thesis, the insights, detailed behaviors, and components of each member of the diagram will be thoroughly described. ASVv Guidance Navigation Control Water Quality Sensors Raw positions, etc Instructions from Coordinator Motor Velocities WQP Measurement Information Movement Feedback from Real World Figure 3.3 Simplified on-board systems of each ASV. In this thesis, these systems are considered to be available and fully functional. 3.1.4 Objective Functions Multiple WQPs must be obtained simultaneously, and this mission can be seen as a Multi-Objective Optimization equation, which typically has the form of: optimize { s=S ∑ s=1Fs(x)} s.t. x∈X (3.4) where Fs(x) is one of the S functions to be optimized (these functions are not the unknown functions fs(x) ), i.e., F1(x),F2(x),...,FS(x) . Generally, there is no solution that simultaneously optimizes all functions, since the objectives are counterbalanced, that is, the value that optimizes Fs(x)is not equal to the value for optimizing any other function Fs(x). Each function Fs(x) , rather than each WQP function, corresponds to a summarized error value of the surrogate model obtained of a real WQP s in a particular water body. This thesis considers each function Fs(x) as a function that evaluates the average difference 28 Chapter 3. Methodology Centralized Coordinator Server ASV1ASV2ASVV Locations where to measure Measurement Values Surrogate Water Quality Parameter Models Figure 3.4 System Overview comprehending a Centralized Coordinator, ASVs, which in turn has GNC and WQP Sensor systems, and the expected output (reliable surrogate WQP Models). between a real or objective function fs(x) and the surrogate or approximate regression ˆ fs(x), and this value is considered a performance metric. Performance Metrics To estimate the model error, the mean squared error (MSE) is usually taken as the loss function. MSE measures the quality of a surrogate model w.r.t. the real model in terms of squared differences between them. MSE yields an absolute error value, which cannot be used directly to evaluate multiple MSEs simultaneously because, as shown in Table 3.1, the real values are not normalized and therefore an average value is not directly available. Consequently, another metric will be used to account for this normalization. The Coefficient of Determination is a regression analysis estimator that evaluates how well independent variables can explain the dependent variable. Moreover, this measure, often named R2 or R-squared, also depends on the square difference, as does MSE, and additionally depends on the variance of the data, which acts as a normalization factor. Then, a single average of R2 values can be obtained for multiple different WQP models. R2has the form of: R2 s(fs(x),ˆ fs(x)) = 1−∑(fs(xi)−ˆ fs(xi))2 ∑(fs(xi)−¯ fs(x))2(3.5) In this equation, ¯ fs(x) corresponds to the mean value of the real data. Since the equation evaluates the proportion between the sum of squared errors (i.e., the residual sum of squares) and the total sum of squares, the value can vary between (−∞,1.0] , because the proportion, which acts as a variation of the MSE (which improves towards zero), is subtracted from the value of one. The system will be useful when the average R2 of multiple surrogate models is closer to one. 3.2 Bayesian Optimization Framework for Monitoring 29 Multi Objective Optimization Since R2 is the selected performance metric, a maximization is obtained for multiple target functions, which implies that the optimization procedure (Eq. (3.4) ) is to obtain the highest values for each R2 . Furthermore, the maximization of each Fs(x) = R2 s(fs(x),ˆ fs(x)) is not done directly, but by fitting the surrogate models according to the measurements that were recently evaluated. In that sense and with regard to multi-objective optimization, after fitting the surrogate models, a linear scalarizing a priori method is applied, with equal weights, so that the solution corresponds to an average of the R2 values. Definitely, the optimization procedure, obtaining better surrogate models, is done indirectly (through the proposed monitoring system), whereas the performance evaluation is done using the obtained surrogate models and a comparison of them with the real models. 3.2 Bayesian Optimization Framework for Monitoring With the environment appropriately defined both mathematically and according to the water monitoring application, the next step is to formally define the system, which applies the Bayesian framework to obtain measurement locations. For simplicity and without loss of generality, the present and future sections will consider only one WQP and one ASV, such that a bottom-up architecture can be appropriately defined afterward. Originally, the idea of BO focuses on optimizing (finding maximum/minimum values) an unknown function, which is not only costly to evaluate, but its first and second derivatives are also unknown to the decision maker and the system [ 79 , 65 ]. This optimization is done sequentially, through the usage of utility functions based not on the real unknown function but on a surrogate model that yields an approximate function, so that in real terms the optimization is done to the surrogate model instead of the real one, while the surrogate model is sequentially expected to resemble the real model. Hance, BO has the general form of: x∗=arg min x∈X f(x)(3.6) which implies that the optimization goal is to find the minimum location for the black-box function, i.e., x∗is the measurement location that optimizes the WQP function f(x). Considering a proposed measurement location x∗ and the fact that f(x) is unknown, the real function is replaced by the regression produced by the surrogate model ˆ f(x) . Then, to update the belief, the function is evaluated (WQP is measured) at this location x∗ . Now, the system can incorporate the newly obtained information, update the surrogate model, and, accordingly, propose a new location for measurement x∗ . This is why BO is a sequential method and, since it depends on surrogate models, it also uses utility functions that may depend on stochasticity. The right-hand side of Eq. (3.6) can be replaced with these components, considering that the utility function α(·) depends directly on the surrogate model. 30 Chapter 3. Methodology 3.2.1 Surrogate Model The BO paradigm is to define a location that optimizes an unknown function under a lack of information or uncertain beliefs. For that matter, surrogate models are used to infer the behaviors of the true black-box functions. These models can be both deterministic and stochastic. The latter is the most widely adopted because uncertainties can be used efficiently to balance exploration/exploitation through specific utility functions. Stochastic models yield not only the predicted mean value but also the related uncertainty. In the literature, ˆ f(x) corresponds to the mean value of the stochastic model but is more frequently referred to as µ(x) , while uncertainty is referred to as σ(x) . Gaussian Processes are often used as probabilistic surrogate models, since they include prior beliefs (such as expected behavior and covariance) and fit the regression through posterior updates (using evidence). Moreover, GPs are relatively easy to understand because of their derivation from normal distributions. GPs are formally described in the next section (Section 3.3). 3.2.2 Utility Function Utility functions often use the GP regression to define optimal evaluation locations. Generally, these functions consider both the mean and the uncertainty to acquire a new location to measure. It is evident that this proposed location is x∗ , and that the complete implementation of the BO procedure is redefined from Eq. (3.6) as follows. x∗=arg max x∈X α(µ(x),σ(x)) (3.7) considering that the mean and uncertainty are obtained using a fitted GP. Moreover, the procedure is now to maximize the function, mainly because, generally speaking, these functions are designed that way. Note that if we consider a weighting factor for both arguments of the function α(·) , the acquired measured location x∗ can be based exclusively on the mean or uncertainty. When considering only the uncertainty, the procedure is often named Active Learning, in which the idea is to measure the most unknown location to learn the model as accurately as possible [ 81 ]. In movement-free environments, this procedure would be preferred, since the mission is to obtain fair approximations of the unknown environment. However, since vehicles must travel from location to location, the system can also make efficient use of the mean value to obtain the measurement locations. Generally, the specific function is named Acquisition Function. An extensive study on Acquisition Functions can be found in Section 3.4. 3.3 Gaussian Processes as Surrogate Models Gaussian Processes are a generalization of the Gaussian probability distribution [ 20 ]. Based on multivariate Gaussian distributions, GPs can fit data with ease (some examples include [ 19 , 67 , 12 ]) through the definition of its 2 main components: the mean function µ(x) and the covariance function k(x,x0) , both of which are expectations of values according to a 3.3 Gaussian Processes as Surrogate Models 31 function, more formally: µ(x) = E[f(x)] k(x,x0) = E[(f(x)−µ(x))( f(x0)−µ(x0))] (3.8) Note that these functions are defined according to the real process f(x) , but since the process is unknown and its evaluations can be noisy, it is safe (and almost necessary) to define the target function also as stochastic, with this an approximation of f(x) can be obtained with the defined mean and covariance functions, according to [20]: f(x)∼GP(µ(x),k(x,x0)) (3.9) The equation above implies that ˆ f(x) = GP(µ(x),k(x,x0)) , so the regression of the surrogate model in this definition yields a stochastic-dependent variable. In conclusion, GPs can be described as a way to obtain a Gaussian distribution for each input variable. In the case of this thesis, this means that a GP can obtain the expected value of a WQP (and its standard deviation) for each location inside the water body. One of the main advantages of GPs is that they support the fit of the model through modifications (marginalization and conditioning) of their parameters that will maintain their closed behavior, i.e., the result of applying these particular modifications yields another GP [ 82 ]. Since the confidence interval encloses a volume, the 2D GP regression is often represented by two 2D heat maps showing the mean and standard deviation functions of the model. The dependence between the input locations is modeled by the covariance function, which will be defined in the next subsection. Subsequently, the process of obtaining a regression is described in a suitable way. 3.3.1 Covariance Functions Covariance functions, or kernels, are the functions typically used for pattern analysis that are designed to provide a measure of similarity between two random input variables. The relationships values can be grouped into three different ones, but but since the desired GP regression should be able to describe a 2-D map of a supposedly continuous spread of a water quality measure, the studied kernels should be as smooth as possible and present positively correlated values for similar inputs, i.e., close-by locations present similar values of WQPs. Kernels describe every possible covariance between the inputs, and, usually, to ensure handling, this description is modeled as a matrix K , whose elements Ki,j=k(xi,xj) , this matrix is called a Covariance, Gram or Kernel Matrix. For GP kernels, this matrix is always positive semidefinite, which implies a closed behavior, always forming a square matrix. Therefore, to obtain an inferred distribution, the kernel matrix can be augmented to include covariance values between the test locations x∗ with themselves and the covariances of them with the known data input xusing: K=K K∗ KT ∗K∗∗=k(x,x)k(x,x∗) k(x∗,x)k(x∗,x∗)(3.10) 38 Chapter 3. Methodology will be based not only on the GP but also on the current state of the ASV. Acquisition functions are used to determine the projected utility of measuring a specific point within the search space. Usually in BO, an agent seeks to explore during the beginning of the mission and exploit the maximum/minimum values towards the end. There are several widely used acquisition functions, each balances exploration/exploitation according to the design of the function. This means that the balance of exploration/exploitation is done automatically with the usage of the functions. However, most acquisition functions include an explicit exploration/exploitation parameter ξ , which manually weights the parameter, which in our case should always favor clearance of uncertainty over optimization. All selected acquisition functions will undergo test scenarios and the one that selects the best measurement locations on average will be selected to be part of the final system. Before defining some acquisition functions, it is better to define the normalized optimum response of the surrogate model, which is named Z . According to [ 79 ], with the current minimum (best) value measured f(xB),Zhas the form of Z=(f(xB)−µ(x)−ξ σ(x)if σ(x)>0 0otherwise (3.26) which implies that Z will have higher (positive) values for high negative mean values or high covariance values. Having Z defined, some of the most well-known acquisition functions in a BO setting are described below. •Probability of Improvement (PI): It computes the probability of obtaining a value that is better than the current best minimum value. It is based on the Cumulative Distribution Function (CDF) Φ(Z) of the normalized surrogate model Z . The probability of improvement of a location can be obtained simply using the following. PI(x) = Φ(Z)(3.27) Since Z is large for an optimum minimum point, it is expected that PI will always prefer to exploit the black-box function once the system evaluates it near a local/- global optimal region, because the CDF is proportional to Z . This statement is true for several reasons, in the context of this thesis, including having a limited bounded box region, expecting a normal distribution that will limit the value of µ(x) that we can evaluate, and the fact that WQPs are expected to have simple, continuous, smooth behaviors, despite being multimodal. In that sense, PI is usually better for unimodal behaviors (or for very few optimum locations). •Expected Improvement (EI): It improves the PI function using predictive uncertainty to prevent exploitation bias that occurs when only the probability of improvement is considered. Therefore, EI uses not only the CDF but also the Probabilistic Density Function (PDF). EI returns a measure of optimization beliefs based on the 3.4 Acquisition Functions 39 current CDF Φ(x)and PDF φ(x)of the standard normal Z. EI(x) =      0if σ(x) = 0 (f(xB)−µ(x)−ξ)Φ(Z) +φ(Z)σ(x)else (3.28) The exploration/exploitation balance is implicitly found in this AF. The predictive mean (found in the first term of 3.28) manages the exploitation weights, while the predictive uncertainty σ(x)provides more exploration weight as it increases. •Scaled Expected Improvement (SEI): It is an extension of the EI method proposed in [ 85 ]. In this function, the EI (x) is scaled to favor the selection of points where the improvement has a small variance, but the predicted value is expected to be high. The scaling factor V[x]is calculated as V[x]2=k(x,x)[(Z2+1)Φ(Z)+Zφ(Z)]−EI(x)2(3.29) The SEI is expressed as SEI(x) = EI(x) V[x](3.30) Active learning should also be introduced in this subsection since it is a straightforward and simpler approach to obtain unknown models based on exploration. But it will be seen in the next sections and chapters that it is not guaranteed to provide the best results in a constrained movement scenario like monitoring through ASVs. Active Learning considers Eq. (3.7) weighting a complete bias towards the uncertainty of the model, such that the next optimal location corresponds to the location in which the uncertainty is the maximum, or, more formally: x∗=arg max x∈X σ(x)(3.31) which can be considered as one of the possible AF to be used by the system. Note that the usage of this Acquisition Function (AF), as well as any other, will be selected according to the constraints of the ASVs, which will be discussed in later sections. Now, with all the AFs presented, it is important to state that the general idea is to select one of the functions according to their expected utility; i.e., exploration-based functions should be used when the uncertainty of the GP needs to be minimized, while exploitation parameters/functions should be favored whenever a certain value of the function is being searched. These AFs provide useful results in BO when used to optimize costly functions. However, they present some limitations related to their application to the monitoring task of an ASV. Therefore, additional work is required to adapt them so that BO can be applied in large-scale water environmental monitoring tasks. The process of obtaining a new measurement position in this thesis differs from common approaches, since the assumptions of the search space are also different. Water bodies are generally non-convex sets, and the definition of the space is also limited to where the water 40 Chapter 3. Methodology or the terrain is present. This constraint leads to the recognition that the monitoring space is not continuous on the shores of the lake as water quality information cannot be obtained in the terrain. With these assumptions, commonly-used fast optimizers like the l-BFGS or the Newton-Conjugate Gradient cannot be used, as they only behave correctly whenever the search space is continuous. l-BFGS-B works great for general BO approaches; however, in the general scenario, a box constraint is present and every point inside of the boundary is available for evaluation. This, of course, is not the case in this thesis, so other methods must be addressed to obtain the best possible next measurement location. 3.4.1 Proposed Acquisition Functions To overcome the limitations of classical AF for the monitoring task of an ASV, new adaptations are proposed to include mechanisms that encourage the selection of points that are close to the current position of each vehicle. Three different adaptations have been designed and tested; each focuses on achieving optimal measurement locations near the current position of the ASV. The main idea is to take advantage of the long distances traveled from one measurement location to another. Consider the following figures as shown in Figure 3.5 in which there exists a previous measurement already done at locations marked as yellow triangles, the last location that coincides with the current location of the vehicle marked as a red triangle and the measurement locations marked as blue crosses according to one of the Adapted Acquisition Functions that are described as follows. Adaptation functions •Split Path: During long distances, the ASV can take several samples before updating the surrogate model. This method calculates the new target location according to the classical AFs and the current ASV position. Then, a number of measurements is taken within the segment from the vehicle’s position and the target point. These locations are equally distanced between each of them by a constant l , which defines the average distance that the ASV should travel to perform a new measurement and is different from the hyperparameter length scale ` . Figure 3.5a shows the waypoints defined by using AF with the split path adaptation. The vehicle should visit all goals and perform measurements on each point, but will only fit the surrogate model and calculate the next optimal location whenever measurements on all points defined by this method have been performed. The real distance traveled between the waypoints will generally be equal to l , except on the last path or whenever the straight line between locations crosses obstacles. •Truncated Path: In this method, the direct path to the next target location according to the AF is “truncated” after a distance l . The idea is to take into account the mobility restrictions of ASV, so it should not travel long distances before updating the surrogate model. Notice that the next optimal measurement location is the same as the first in the split-path method, since it is defined by the above-mentioned AF. In fact, Figure 3.5b shows a measurement location that is the same as the first in the previous subfigure. However, in this adaptation the distance traveled by the ASV is truncated (hence the name), reaching a suboptimal location. Whenever the ASV reaches the computed location and updates the data, a new goal location is obtained. 3.4 Acquisition Functions 41 This new location is generally not the same as the second measurement location using the split-path adaptation due to the new information provided. •Masked Path: This method weights nearby locations using a simple gradient according to (3.32), where x corresponds to the locations on the search space and pv the position of the ASV. This process is carried out by multiplying the mask by the result of the classical AF, as shown in (3.32). This method allows for exploiting the current position of the ASV in order to obtain more information with less movement. M=exp −(||x−pv||) l(3.32) Due to the exponent, the rapid decay that this mask applies will give greater values to locations that are closer to the current location of the vehicle. Figure 3.5c shows the mask that is applied to the AF, and the resulting AF is shown in Figure 3.5d. A summary of the proposed adaptations is included in Table 3.3, showing that they depend only on a tunable hyperparameter l . It is important to note that each has a distinctive characteristic with respect to updating the surrogate model. The split path-adaptation will perform at least one or more measurements before fitting the Gaussian model. The two others will strictly fit the model after each measurement. The truncated-adaptation provides the information to the ASV in its path to the maximum, so it can quickly update the model. Finally, masked-adaptation favors the maximum of nearby locations. Table 3.3 Summary of properties of the proposed adaptations.. Adaptation Distance between new Measurements done data locations and ASV before GP Fitting Split path l,2l,...,nl,||pv,pm|| n+1 Truncated l1 Masked undetermined 1 Definition of maximum distance for adapted measuring In this thesis, it is preferable that the length of the new segment l is dynamic and related to the length scale ` of the fitted (posterior) GP ( l∝` ). This rationale is based on the premise that the length-scale ` defines the maximum distance between the inputs ( x∈X ) that allows these inputs to influence each other. Therefore, there exists a dependence between the optimal distance between measurement locations and the underlying objective function, i.e., large distances may fail to model the objective function, while very short distances can provide redundant information. Hence, since the length-scale is the only parameter that is updated taking into account the real evaluations to fit the model, an optimal distance between measurement locations can be obtained using the length-scale. The use of length scales to limit the distance traveled is suitable and interesting because, for example, in the SE kernel, ` captures information about the spatial frequency response. 42 Chapter 3. Methodology (a) (b) (c) (d) Figure 3.5 Subfigures (a), (b) and (d) correspond to examples of different results using the 3 different adaptations proposed in this thesis for the same initial conditions. Sufigure (c) is the mask applied in the masked adaptation. The images corresponds to a work previously made by the author. To prevent loss of information, the parameter λ is set to further limit the traveling distance. Finally, the value of lis always given by the following equation. l=λ×`(3.33) 3.5 Multi-Water Quality Parameter Estimation At this point, a single vehicle could efficiently obtain a WQP map model just by using an appropriate kernel, and a good acquisition function that is adapted to the needs, but since there exist multiple WQPs, the next step consists of selecting a measurement location 3.5 Multi-Water Quality Parameter Estimation 43 according to multiple objectives, i.e., solving a Multiple-Objective Problem. The objectives are defined following the aforementioned BO approach. 3.5.1 Multi-objective problem Multi-Water Quality Parameter Estimation is done by slightly modifying Eq. (3.6) , and consequently Eq. (3.7) , to account for multiple fs(x) . Evidently, this implies that there could be different GP kernels that would fit the different WQPs since their behaviors are not necessarily the same and do not have the same scale factor. Furthermore, the acquisition functions could not necessarily be the same, implying different αs(·) functions, considering the different WQPs. Nonetheless, since only one ASV is available (in this section of the chapter), only one measurement location should be selected considering this MOO setting. Regarding the adaptation, is is not considered as an objective, but as a posterior modification of the selected measurement location, therefore it is not included and the acquisition functions are the classical without modifications. The generalization approach consists of determining an appropriate location x∗ so that a set of surrogate models can be efficiently approximated simultaneously. Next, two ways for solving the multi-objective problem are presented. On the one hand, Multi-Function Estimation (MFE), which linearly combines every AF. On the other hand, the second method finds solutions according to a Pareto-based Multi-Objective Optimization approach. 1. Multi-Function Estimation (MFE): For simplicity, the MFE AF can be defined as a composite of different AFs, or a general AF that takes into account the different functions, such as the work in [ 71 ]. The proposed approach consists of composing different AFs through a fusion procedure. Based on [ 71 ], two different AF Fusions are defined, i) the decoupled and ii) the coupled methods. The only requirement for these methods is that the acquisition functions are normalized, which is normal if the functions utilize the normalized optimum response (Z). a) Decoupled Evaluation: It is designed to optimize one single WQP map model at a time, so that different objectives are optimized in different steps. The expression responds to select the next measurement location x∗ as the argument of the maximum of the different AFs (eq. 3.34). For example, if all AFs weight the different locations of the set X according to the uncertainties σs(x) of the different models, the measurement location will correspond to the position where one of the AFs has the highest uncertainty value. The decoupled method expression is shown below: x∗=argmaxx{αs(x)},s= [1,2, ... , S](3.34) b) Coupled Evaluation: It is a fusion that acknowledges the importance values of all AFs. Therefore, AFs contribute equally to a general AF. It consists of adding all the AFs. Thus, the next measurement location x∗ is calculated as the argument of the maximum of the sum of the AFs, as shown in eq. (3.35). This location x∗ 44 Chapter 3. Methodology will be selected as the location where the different AFs will have their values maximized as a combined set. x∗=argmaxx{ S ∑ s=1 αs(x)}(3.35) 2. Multi-Objective Optimization (MOO): The AFs can also be considered objectives under a MOO framework. Therefore, the optimality of a solution is based on the Pareto dominance over the objective space. A solution, which in our case is a candidate measurement location, is said to be a non-dominated solution when no other solution can maximize one of the objectives without affecting (degrading) the others. The mentioned solution is called Pareto efficient, and the set of Pareto points is called Pareto (efficient) Set. Therefore, the proposal is to use a Pareto-based approach to obtain solutions (measurement locations) that balance the multiple objective functions (acquisition functions) simultaneously. A Pareto set is obtained using the Non-dominated Sorting Genetic Algorithm II (NSGA-II) [ 86 ] as provided by the pymoo python package [ 87 ]. Being a Genetic Algorithm (GA), this method requires the usual GA definitions: individual, population size, number of generations, and the genetic operators: crossover, mutation, and selection functions. • Individual representation: Consists of a single 2D point, which indicates the location of a measurement location. • Fitness function: Consists of a vector of AF values at the individual measurement location, such that every value is directly obtained from the selected AF α(·)considering the different models to be obtained. • Genetic Operator, Mutation: Polynomial Mutation [ 87 ] has been selected so that individuals mutate locations according to a probability distribution, effectively generating offspring that explore their surroundings to find local optimal values. • Genetic Operator, Crossover: Simulated Binary [ 87 ] is selected to exchange information between individuals. It works in a similar fashion to that of the Polynomial Mutation operator, creating a crossover that takes into account the data of two individuals and combines them. • Genetic Operator, Selection: Because of NSGA2[ 86 ], Tournament Selection [ 87 ] has been used due to its simplicity and its usefulness for faster convergence. This operator selects the best individual out of a group according to the fitness function. Since the fitness function is multi-objective, the complete process is based on Pareto efficiency. The result of this method is a set of optimal (non-dominated) measurement locations. Finally, since the method provides a set of possible measurement locations, the closest non-dominated solution to the vehicle is selected as the next optimal measurement location x∗ . This selection is done based on the assumption that measuring on any of 3.6 Multiple Autonomous Surface Vehicles 45 the Pareto efficient locations will yield the same performance. Therefore, selecting the closest location is a suitable approach considering that the measurement should also be energy efficient. 3.5.2 Multi-function adaptation Whenever a location is selected to be the optimal solution for a future measurement, the next step is to adapt the location so that it can still account for energy consumption; this is to apply one of the adaptations found in the previous section. Subsequently, the length l is changed so that it corresponds to one of the different length scales, more formally: l=λ×min{`s}(3.36) 3.6 Multiple Autonomous Surface Vehicles This thesis proposes a system that monitors multiple WQPs using multiple ASVs. Since the definition of monitoring multiple WQPs has already been described in previous sections, the usage of multiple ASVs is left to be defined. The idea is that each ASV cooperates with the mission according to the collective beliefs and their locations in selecting the measurement location. Cooperation, on a broader scale, can be divided into two categories: i) intentional and ii) unintentional cooperation strategies. Each of them presents different advantages and may be preferable depending on the environment in which the agents operate. In an active or intentional cooperation setting, the agents share their current goals and information, decide on steps such that every other agent will obtain benefits from the actions and there is a constant interchange of information so that there is a global balance in energy usage/mission accomplishment. It is generally preferred when agents need to exchange a lot of data, cannot fully operate without constant communication, and there is a need for balance in terms of energy usage. On the other hand, unintentional cooperation describes agents that collectively pursue a mission with actions that are not deliberately joint or on purpose. In the case scenario of this thesis, this strategy would force ASVs to work on their own with the objective of completing the mission without an intended joint effort with other agents. Note that unintentional is not the same as unintended; the last adjective describes that the cooperation was not planned, but resulted as an outcome. In this thesis, ASVs are expected to cooperate, but their actions are not affected or do not consider the actions of other vehicles. The main advantage of this strategy is that an asynchronous execution can be planned and that it does not matter whether another vehicle completes its mission or is added to the fleet after another agent has already defined a measurement location goal. Despite the chosen cooperation strategy, the plan is to obtain as many WQP measurements as possible using the available ASVs. In this sense, a centralized coordinator is needed to collect the data gathered by the vehicles. The ASVs will work (move and perform measurements) independently, they will communicate to the coordinator their locations, and the values of the WQPs according to their onboard sensors. This implies that the ASVs would only need to receive mission measurement locations instead of receiving 46 Chapter 3. Methodology all the data from other vehicles, computing a local Gaussian process, applying the BO equation with an AF, and obtaining a measurement location. In a sense, each ASV is unaware of the BO process, only performing Path Planning from one location to another, performing the measurement, and communicating the results afterward. The unintentional straightforward cooperation strategy is to relegate the decision making solely to the following: i) the centralized coordinator, ii) the available WQP data at measurement location generation time, and iii) the current locations of the different ASVs. On the other hand, intentional cooperation will also consider the current missions of other ASVs (the location where they are heading) and the data that they would gather. The unintentional method will definitely work better with fewer agents available that are capable of joining or leaving the mission at any time, while the intentional method will work better for missions where there is certainty that the starting ASVs will endure throughout the whole mission, while the addition of ASVs in the middle of the mission would be counterproductive in terms of redundancy and energy efficiency. Another important factor is initialization; normally, ASVs would start from the same harbor, but this is not necessarily the best approach, since vehicles in the same location would be too close to each other to provide useful data (redundant measures). A way to overcome this issue is to efficiently distribute the initial vehicles on the surface of the water body. 3.6.1 Initial Vehicle Positioning A mathematical approach to solve this problem is to distribute the vehicles in the perimeter of a radius circle r centered at the center of the water body, evenly separated angle-wise, in an attempt to have each vehicle aimed at covering different regions. However, for real scenarios, it is not desirable that the ASVs start the mission in the center of the water body. Therefore, the initial positions are relocated to the shore of the water body taking into account the spacing between vehicles, having all vehicles evenly separated angle-wise but with different distances with respect to the center of the lake according to the shoreline. 3.6.2 Region Partitioning using Voronoi Diagrams Having the vehicles at “optimal” initial locations, the next step is to partition the feasible space according to the number of available ASVs. In this thesis, Voronoi Diagrams (VDs) can be effectively used to do so and obtain evenly distributed regions. VD is a partition model of a n -dimensional space into regions of the same dimension according to a distance rule. In this work, the 2D space of the map M is divided into v regions, that is, one region for each of the V vehicles available. Each region has a generator, which corresponds to the position or location pvof the vth ASV. The VD expression has the form Rv={x∈R2|d(x,pv)≤d(x,pw)∀v6=w}(3.37) where d(·) is assumed to be the Euclidean distance function. Using VD as a partition system is efficient and appropriate, but have some drawbacks that need to be addressed. Two generators cannot be at the same location. This is completely avoided if the initial locations are not close to each other. The outermost generators always have an infinite surface area, because the distance to a location x will be smaller than that of any other 3.7 Proposed Multi-ASV Multi WQP monitoring system 47 generator. To overcome this problem, it is important to consider only the feasible locations in M . Finally, these regions are updated as needed, which means that every time an ASV needs a new measurement location, a new VD is computed, the new VD will also comply with the former definitions so that it is always valid. The use of VD not only helps to appropriately and efficiently divide the regions according to the total number of vehicles available, but also strongly contributes to scalability and robustness. The periodic recalculation of measurement VD regions for each ASV provides an easy and reliable mechanism for robustness, because ASVs could be included and removed from the system at any time. Thus, the location of the measurement will be adapted to the number of available ASVs. 3.7 Proposed Multi-ASV Multi WQP monitoring system In this section, the proposed approach is adequately defined. Figure 3.6 shows the proposed system, where the CC contains a database, as well as three components that work together to obtain measurement locations. For the definition of the waypoint or measurement location, the CC system performs three main tasks: i) region generation, ii) surrogate modeling acquisition, and iii) evaluation of acquisition functions. The first step finds the measurement region of the v th vehicle to nearby locations, so that a vehicle cannot find an optimal measurement location outside of its region. This search begins in Step ii), modeling the WQP maps with surrogate models. These surrogate models need to consider any previously acquired data (from all ASVs) and the inferred knowledge (prior models) of the environment, and iii) with the search region defined and the surrogate models fitted, a new measurement location is obtained according to the acquisition functions that should consider the constraints of the ASVs. The three-step process is repeated each time a vehicle needs a measurement location and does not depend on volatile data; therefore, this process can be accomplished asynchronously among vehicles. It is important to note that this behavior helps to alleviate the workload of the CC, as it does not need to continuously transfer data to the available ASVs. In the following subsections, we appropriately define the tasks mentioned. The general procedure of our centralized monitoring system is shown in Algorithm 1. Whenever an ASV needs a new goal (line 4), a Voronoi region is defined using the positions of the available ASVs (lines 5-6). Next, with the available data, the surrogate models are obtained (line 7) and used to find the optimal measurement location within the Voronoi region of each ASV (line 8). Note that this line is generalized with respect to the acquisition function, and that where the subscript s appears, it considers all WQP sensors. It can use MOO using GA or MFE procedures. Recall that parameter l is obtained using Eq. 3.33. Finally, the CC stores new data Di+1←(si+1,pi+1)=(sv,pv) when received/available and proceeds to fit the models (lines 10-11). 54 Chapter 4. Implementation WQPs. Moreover, specifically for WQPs, the benchmark functions should be multimodal, continuous, multidimensional, and deterministic. Finally, real WQPs are not expected to behave always like the mentioned benchmarks, and it was more useful to create custom benchmark “ground truths”. This is because a WQP is expected to have many peak locations, to be continuous everywhere (except on the shore), to depend on at least two independent variables (x,y) , and to have an underlying ground truth value to which it could be compared. The Shekel function (SF) is a rather unknown benchmark that can be adjusted to provide maximum locations and has gentle gradients. Additionally, the SF can have any number of maximum values, which fairly simplifies the setup procedure to obtain the ground truth for each simulation run. SF has the form of: fShekel(x) = M ∑ i=1 1 ci+∑N j=1(xj−ai,j)N(4.1) where ai,j and ci are the elements of two sets of matrices A and c . The matrix A has the size of M×N , where M is the number of maximum locations, which are N dimensional. The matrix c is a vector M×1 that defines the inverse importance value of the maximum locations. The SF provides a parametric function capable of defining multiple optimum locations in a N-dimensional space. Different ground-truth maps can be generated using Eq. (4.1), each can represent different maps of water quality parameters at different times. They can be post-processed so that they can fit the map M and µs(x) = 0 , σs(x) = 1, making the GP fit easier and faster. This post-processing corresponds to subtracting the real obtained values with the expected mean of the measured WQP and dividing by a factor equal to the standard deviation of the measured WQP, effectively normalizing the map, assuming the Central Limit Theorem. Examples of the generated ground truths can be seen in Figure 4.3. In either way of choosing WQP models, the different simulated agents should travel from point to point until their battery is fully used, and they performed multiple WQP measurements. 4.3 Autonomous Surface Vehicles Regarding ASVs, the implementation was done so that the use of simulated vehicles and real vehicles is not affected by the system. In that sense, the same components of high-level systems have been designed and implemented for both real and simulated ASVs. In what follows, the descriptions of the components will be detailed for the real ASVs, but they also work for the simulated ASVs, unless explicitly described, which refers only to one of them. 4.3.1 Guidance, Navigation & Control (GNC) The Guidance Component is designed to apply Eq. (3.7) , under the assumption that a single next measurement location pm+1=x∗ is required for a single vehicle Vv that needs a target location. Evidently, depending on the number of agents present at that time, the available locations x varied according to the region obtained for that vehicle in that scenario Rv (a complete region if only one vehicle is present) and the acquisition function α(·) 4.3 Autonomous Surface Vehicles 55 0 200 400 600 800 0 200 400 600 800 1000 1200 1400 0 200 400 600 800 0 200 400 600 800 1000 1200 1400 0 200 400 600 800 0 200 400 600 800 1000 1200 1400 0 200 400 600 800 0 200 400 600 800 1000 1200 1400 Figure 4.3 Four examples of Shekel functions with random number of maxima points at random locations scaled to a map model of the Ypacarai Lake. selected. In the case of simulations, multiple α(·) were tested, so this can vary whether it considers Expected Improvement, Probability of Improvement, etc. for one, the WQPs or AF fusion with multifunction estimation or multiobjective optimization for multiple WQPs. After a measurement location has been obtained and sent to an ASV, the global and local path planner components consider the current location and the goal of obtaining an obstacle-free path; this is done using RRT*. These obstacle-free paths were obtained for the map model Minstead of the real world map, as they were easier and faster to solve. During the mission, the navigation system of an ASV provided location data to the CC as well as other components of the ASV, so efficient paths were obtained and efficient path following could be executed. As mentioned above, this component was executed by an embedded autopilot system such as Navio2. This autopilot system also executed the control system. 56 Chapter 4. Implementation 4.3.2 ASV Constraints Regarding the main ASV constraint, battery consumption, for testing purposes, the maximum distance was drawn from the theoretical maximum battery available according to the tests performed in [ 91 ] with the average velocity considered in the same work. As a result, the simulated ASVs obtained the best possible maps. In real experiments, the battery will be used until a certain threshold level that will automatically return the vehicle to launch, or RTL as described in the Navio2 autopilot system. 4.3.3 Water Quality Paremeters Sensor System WQP sensor systems are the main part of the ASVs, so they were always available for use. If any WQP system of an ASV is supposedly malfunctioning, the ASV will no longer be considered during the mission, so the ASV can return to the shore, be checked, and reenter the ASV system. This is possible by using the unintentional cooperation described in the previous chapter. For real measurements, a certain number of measurements will be performed at the mentioned location, so that the average values can be obtained. In that sense, the system available in the ASVs for real experiments corresponds to Libelium Smart Water [ 92 ] (Figure 4.4), which has the same physical and chemical sensors as shown in Table 3.1. On the other hand, for simulated ASVs, a measurement value was drawn directly from the simulated map, and a random i.i.d. noise was added to account for noisy sampling. The CC is an online server that acts as both a decision maker and a visualizer. All data that arrive from the vehicles were used to produce output, such as surrogate models, methods for controlling the state and health of the vehicles, etc. Figure 4.4 Water Quality Parameter Sensor system, Libelium Smart Water, connected to the ASV. The sensors are connected to the bottom of the embedded system and communicates to the main on-board system through an USB cable. 4.4 Implementation using Python 57 4.3.4 Centralized Coordinator Server The CC is the only system that worked directly with Bayesian optimization, Gaussian Process, Acquisition Functions, and Voronoi Diagrams. Real vehicles will receive and send data via 4G, assuming that the mobile phone signal is available anywhere on the water body. The data will be sent without processing them, so only a limited amount of processing is done on the onboard computers to save energy. In that sense, the energy usage of vehicles was considered to only power the motors to perform movement. 4.4 Implementation using Python During simulations, the main tool for comparison is the Coefficient of Determination R2 . The idea behind always considering this performance metric is to have a systematic method for comparing not only the component selection but also to observe and analyze the joint operation of multiple components such as Multi-Objective Optimization or Multi- Agent system operation. When a ground truth is available (which only happens during simulations), R2 could always be obtained. Therefore, this decision is plausible and valid. Considering the framework of BO already described, the implementation is to simply code in Python language Eq. (3.7) . Then publicly available libraries are used to implement parts of the system in Python language; in some cases modifications of the original library were done so that the implemented code could work seamlessly within the system. The usage of the main libraries involved in the implementation are mentioned next. •Numpy: https:// numpy.org/ is a mathematical tool used for data handling, all data are mainly handled by this package. All vectors, positions, and most operations are done using this package. Some examples include matrix multiplication, battery usage, and path planning. The libraries described below also use numpy to perform their operations. • Scikit-Learning: https:// scikit-learn.org/ stable/ is the Python package that handles kernels, such as SE, Matern, and others. Additionally, it implements Gaussian Processes, such as data processing, fitting, and even hyperparameter optimization. It is very useful for noise-free sequential decision making. However, for noisy evaluations, the package is not fully implemented. The hyperparameter optimization in this library implements Eq. (3.24) with the observed data and provides the mean and uncertainty of the fitted GP using equations. (3.19) and (3.20) , using Cholesky decomposition where possible. • Scikit-Optimize: The bases of classical acquisition functions are drawn from this package https://scikit-optimize.github.io/stable/. To enhance response time and fully utilize the information given by the AFs, all utilized AFs are reimplemented as new classes with added functionality. • Pymoo: When multiple acquisition functions are to be computed and solved as a multi-objective optimization strategy, the library pymoo, https:// pymoo.org/ , is used. It is a very high-level library that, through a problem definition, could provide a Pareto-efficient set that is supposed to be part of the Pareto frontier. 58 Chapter 4. Implementation • Scipy: https://scipy.org/ is a high-level library for scientific computing. In this thesis, Voronoi diagrams are obtained directly using this library. However, copies of the locations of the vehicles are translated and mirrored to obtain a box partition of the region, since scipy only offers Voronoi computing considering non-restricted regions. • Matplotlib: All figures regarding the results in this thesis were generated using this library https:// matplotlib.org/ . It is a simple library for creating visualizations. In addition, this library is used for visualization of the resulting WQP maps. Since the objectives of this thesis are very specific, additional code was written for the adaptation of the acquisition functions and multi-water quality parameter estimation using coupled and decoupled fusion strategies. Both are available online. These pseudolibraries are largely based on the libraries previously mentioned and work, to the best known knowledge of the author, only with them. In the next chapter, the simulation and experimental results are shown as a process in which, first, the basic components of the system were systematically selected according to their performance metric. 5 Results Taking action today is worth it. Kurzgesagt – In a Nutshell, 2022 In this chapter, the results of the proposed method. based on BO for monitoring a simulated aquatic environment are described. The followed roadmap consisted of evaluating different components, selecting the best method, and building the system from the ground up. Comparisons have also been made with related monitoring systems to show its strengths. 5.1 Simulation Setup This section defines the procedures and parameters for the evaluation of the proposed method. The simulations have been conducted in a 6-core 4.6 GHz processor computer with 16GB RAM. In the simulations, to ensure a fair comparison among different components, the simulations have been performed with the same random seeds. Simulated Water Quality Model using Benchmark Functions The WQP maps shown in Figure 4.3 are examples of ground-truths from the test configuration. The simulated WQP value of each location is taken from the result of a scaled standardized shekel function. For the randomization of the shekel function, random values were defined for the input matrices. More precisely, the number of peaks or maximum values was defined as a random integer between two and ten, the maximum altitude was defined to be a real number of minimum values of five to ten, while the locations of the peaks were distributed uniformly random across the map matrix M. Theoretical maximum energy available According to Morel et al. [ 91 ] and Eq. (3.3) , a maximum travel distance can be obtained considering the average velocity ||~v||=1.6[m s] , the maximum theoretical energy delivered 59 60 Chapter 5. Results by the battery and the consumption rate provided by the same work. These values are also in accordance with the total time that the ASVs are usually available as it was mentioned in the Chapter 4.2 Implementation. As a side note, since the model map consists of a grid map with square of side d≈10[m], the maximum distance in pixels is 1500 pxs. Dmax =Tmaxenergy ||~v|| e=20000 ×1.6 2.1≈15000[m](5.1) 5.2 Kernel Selection The first selection corresponds to the surrogate model. This thesis considers that GPs are the appropriate model, but has no preference regarding the kernel. For the expected quiet waters of the Ypacarai lake [ 3 ], the prior length scale of the selected kernel should be large enough to ensure smoothness, but bounded enough to allow the kernel enough freedom to adapt to varying maximum/minimum data. The covariance functions to be tested depend all on the hyperparameter ` . It was decided that `∼10% is the length of the search space. For other hyperparameters, if any, the most common values [20, 79] are used. Several simulations have been performed for kernel selection, each test consisting of uniformly drawn sets of M∈[15,25,35,50] valid measurement locations. To account for noisy evaluations, a small random value is added to each measurement so that the noise is i.i.d. with ε=0.01 . For each simulation, a kernel is selected and used throughout the experiment, fitting the data and updating their hyperparameters with the equations described in Chapter 3. All kernels are summed with a White Kernel with σ2 o=0.01 to account for noise. The results using the SE, the RQ and the Matérn kernel are found in Table 5.1, Table 5.2 and Table 5.3, respectively. For each map, 10 different simulations have been performed, and the resulting average R2 is presented. The values found in the last column of each table correspond to the average R2 of the five maps. The results marked in bold correspond to the best average value across different kernels. Table 5.1 R2 values obtained for the Squared Exponential (SE) kernel experiments with a prior length scale `=150 with different settings regarding WQP Map model (10 different simulations per map model) and number of measurements M. MWQP Map 1 WQP Map 2 WQP Map 3 WQP Map 4 WQP Map 5 Avg. 10 0.7458 0.8383 0.6765 0.7265 0.8189 0.7612 15 0.8927 0.9229 0.8261 0.8906 0.941 0.8947 25 0.9564 0.9652 0.9352 0.959 0.9811 0.9594 35 0.9652 0.9769 0.9671 0.9561 0.9871 0.9705 50 0.983 0.9935 0.9826 0.9881 0.9952 0.9885 Note that the generated WQP maps have different values, have their peaks in different locations, and are not correlated to each other. Despite this, for each kernel, the R2 score is similar and has less variance as the number of measurement locations increases. Regarding 5.2 Kernel Selection 61 Table 5.2 R2 values obtained for the Rational Quadratic (RQ) kernel experiments with prior hyperparameters `=150 and α=1.0 with different settings regarding WQP Map model (10 different simulations per map model) and number of measurements M. MWQP Map 1 WQP Map 2 WQP Map 3 WQP Map 4 WQP Map 5 Avg. 10 0.7458 0.8383 0.6765 0.7253 0.8159 0.7604 15 0.8925 0.9225 0.8261 0.891 0.9339 0.8932 25 0.9563 0.9652 0.9352 0.9585 0.9811 0.9592 35 0.9651 0.9769 0.9671 0.9562 0.9871 0.9705 50 0.983 0.9935 0.9826 0.988 0.9952 0.9884 Table 5.3 R2 values obtained for the Matérn kernel experiments with prior hyperparameters `=150 and ν=1.5 with different settings regarding WQP Map model (10 different simulations per map model) and number of measurements M. MWQP Map 1 WQP Map 2 WQP Map 3 WQP Map 4 WQP Map 5 Avg. 10 0.6936 0.7925 0.5959 0.707 0.7952 0.7168 15 0.8709 0.9066 0.7505 0.8863 0.9065 0.8642 25 0.9447 0.9592 0.8949 0.9509 0.9619 0.9423 35 0.9527 0.9613 0.9451 0.9468 0.9703 0.9552 50 0.9794 0.9882 0.9709 0.9851 0.9895 0.9826 the different hyperparameters, Figure 5.1 shows the evolution of different hyperparameters for different kernels. Specifically, Figure 5.1a shows the evolution of the length-scale value ` according to the different number of measurements per WQP Map for the SE kernel. Of course, the length scale value varies according to the map, but it can be seen that the value tends to be a constant. In the case of the other kernels, the same behavior is observed. Moreover, hyperparameter α (the scale-mixture), for RQ (Figure 5.1b) kernels, reaches the maximum allowed values in simulations for some WQP maps, which means it tends to ∞ , effectively, transforming the kernels into SE kernels. The application of the parsimony principle [ 83 ] also applies in this case, since a simpler method of modeling should be preferred. The summarized results using the different kernels can be found in Table 5.4. The values in Table 5.4 correspond to the average R2 considering five different WQP maps. The table shows that, for the Ypacarai Lake, using the SE or RQ kernels provides the best results. Consequently, either of these kernels should be indifferently selected. But since RQ is the generalization of SE and the observed hyperparameter value α reached its maximum possible value, SE is the selected kernel. Moreover, RQ is a more complex expression and provides approximately the same results. For this reason, SE is the selected kernel for the proposed BO-based monitoring system, and it will be used as the prior/posterior model for the next simulations. 62 Chapter 5. Results 10 15 25 35 50 Number of measurements 200 250 300 350 Length scale value posterior Map 1 posterior Map 2 posterior Map 3 posterior Map 4 posterior Map 5 (a) 10 15 25 35 50 Number of measurements 0 25000 50000 75000 100000 125000 150000 Scale-mixture value posterior Map 1 posterior Map 2 posterior Map 3 posterior Map 4 posterior Map 5 (b) Figure 5.1 Hyperparamter evolution corresponding to (a) Hyperparameter length-scale ` for the SE kernel, and (b) Hyperaparmeter scale-mixture α for the RQ kernel . Table 5.4 Summarized average R2 values obtained in the simulations using different kernels with prior hyperparameters as shown and different settings regarding WQP Map model (10 different simulations per map) and number of measurements M. MSE (`p=150) RQ `p=150,αp=1.0) Matérn (`p=150,νp=1.5) 10 0.7612 0.7604 0.7168 15 0.8947 0.8932 0.8642 25 0.9594 0.9592 0.9423 35 0.9705 0.9705 0.9552 50 0.9885 0.9884 0.9826 5.3 Acquisition Function Selection 63 Figure 5.2 shows an example of the prediction generated using a SE kernel with the best prior hyperparameters found in Table 5.4. The first subfigure corresponds to the Ground Truth for that experiment, whereas the second and third subfigures show the mean and standard deviation, respectively, using the measurement values obtained at random locations marked with red diamonds. Regardless of the random distribution of points, the same average result is always obtained using this kernel. This is sensible because when using SE, as a covariance function, one expects smooth behaviors and gentle gradients, as opposed to using Matérn, which creates “noisy” covariance between inputs. The hyperparameter length scale was changed to fit the given measurement locations on each scenario. It has been observed that there were different values of length scale when using only 10 locations, but as the number of measurement locations increased, the variance of the final length scales decreased and in most cases was very different from the prior length scale. This implies that the GP can effectively condition a WQP that behaves like the SF. A careful adjustement of the length scale can be done so that the next experiments start with a better understanding of the underlying black box function, but since different models will produce different length scale and since the real behaviors will not necesarrily always be like a SF, the prior length scale is selected as designed first along with the assumed iid noise. With the kernel appropriately selected, the next step is to select an acquisition function for the system. 0 2000 4000 6000 8000 x [m] 0 2000 4000 6000 8000 10000 12000 14000 y [m] Ground Truth Measurements (a) 0 2000 4000 6000 8000 x [m] 0 2000 4000 6000 8000 10000 12000 14000 y [m] GP Predicted ( x ) Measurements 1.5 1.0 0.5 0.0 0.5 1.0 1.5 2.0 (b) 0 2000 4000 6000 8000 x [m] 0 2000 4000 6000 8000 10000 12000 14000 y [m] GP ( x ) Measurements 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 (c) Figure 5.2 Example of posterior GP regression using the Squared Exponential kernel with `=150 . In this case, with the selected random 10 measurement locations, the fitted length scale changed to a value of `=307.346 using the maximization of logarithmic likelihood. 5.3 Acquisition Function Selection In this section, the results are presented to evaluate the different AFs. Since the idea is to obtain the best-adapted acquisition function for the monitoring mission, a full evaluation needs to be performed considering only one vehicle. This implies that for each acquisition 70 Chapter 5. Results 500 0 500 1000 1500 1000 500 0 500 1000 1500 2000 2500 Figure 5.4 Example of bounded VD for the Ypacarai Lake considering 4 (blue diamonds) initial points and their respective mirrored points (red). 0 2000 4000 6000 8000 x [m] 0 2000 4000 6000 8000 10000 12000 14000 y [m] Other vehicles Current vehicle (a) 0 2000 4000 6000 8000 x [m] 0 2000 4000 6000 8000 10000 12000 14000 y [m] Other vehicles Current vehicle (b) 0 2000 4000 6000 8000 x [m] 0 2000 4000 6000 8000 10000 12000 14000 y [m] Other vehicles Current vehicle (c) Figure 5.5 Voronoi Partitioning Region in action for (a) 3, (b) 4 and (c) 5 vehicles. Note the different regions generated automatically updates and distributes regions according to the number of available ASVs. the different combinations of WQP simulated models used, i.e., the different sensors that were used during a mission, where for each simulation model ID, different runs were made according to different initial locations and number of ASVs. Each Sm corresponds 5.5 Multi WQP monitoring using Multi-ASVs 71 to a different WQP model generated using the shekel function with random peaks and locations. Note that there exist 4 groups of number of sensors involved, each one containing 4 different combinations of WQP models so that generalization was effectively tested. Table 5.11 Different groups of sensors involved. Set Sensors Involved 1 [S0,S1] 2 [S2,S3] 3 [S4,S5] 4 [S6,S7] 5 [S0,S1,S2] 6 [S3,S4,S5] 7 [S6,S7,S8] 8 [S0,S4,S9] 9 [S0,S1,S6,S7] 10 [S2,S3,S4,S5] 11 [S2,S4,S6,S8] 12 [S6,S7,S8,S9] The number of vehicles available were 2, 3 or 4 that had initial locations on the shore of the water body and had a full battery level. The CC used MFE (coupled or decoupled) or MOO (through the Genetic Algorithm) to obtain the best locations for measurements, which were adapted to limit traveling distance. The results are grouped by these strategies. For each of the ASV fleets (2, 3 or 4 ASVs), the mean R2 scores for 120 different simulations using different sets of WQP simulated models are shown in Table 5.12, Table 5.13, for the coupled and decoupled variations, respectively. Regarding this approach, Multi-Function Estimation, no additional hyperparameters need to be determined. Regarding the MOO using GA, several hyperparameters needs to be described. Table 5.14 shows the hyperparameters selected. The code was fully implemented using the library pymoo [ 87 ], and the respective parameters were selected to efficiently explore the feasible search space. In this implementation, the constraints functions are defined according to the map model and the Voronoi region for each vehicle. Table 5.15 shows the results for 120 simulations considering a MOO framework. The table is distributed in the same way as the MFE variations found above. The main difference using this strategy is that a population of candidate measurement locations is evolved to obtain an optimal Pareto efficient set instead of obtaining the location according to an equally weighted approach (MFE coupled) or to an unilateral weighting approach (MFE decoupled). Figure 5.6 shows an example of the MOO in action where there exists 4 ASVs measuring 4 different WQPs simultaneously. In the example, the WQPs correspond to Figure 4.3 found in Section 4.2.1. All Pareto Efficient points that were obtained are within the Voronoi region for the green ASV generated with the current locations of the ASVs, not the last measurement locations. It is remarkable how the proposed approach can efficiently 72 Chapter 5. Results Table 5.12 Results for MFE strategy considering coupled method for different number of ASVs and different number of sensors available. # ASVs # SAvg. R2Avg. M 2 20.928 ±0.069 27.498 ±3.294 30.95 ±0.055 28.673 ±3.26 40.949 ±0.035 28.645 ±3.462 3 20.974 ±0.017 38.922 ±3.546 30.981 ±0.009 40.547 ±3.59 40.979 ±0.013 41.474 ±4.78 4 20.985 ±0.01 51.923 ±3.783 30.989 ±0.011 53.872 ±3.696 40.99 ±0.006 54.096 ±4.701 Average 0.969 ±0.039 40.628 ±10.948 Table 5.13 Results for MFE strategy considering decoupled method for different number of ASVs and different number of sensors available. # ASVs # SAvg. R2Avg. M 2 20.931 ±0.067 27.623 ±3.646 30.943 ±0.057 27.999 ±3.563 40.936 ±0.062 28.174 ±4.312 3 20.975 ±0.016 38.372 ±3.089 30.983 ±0.01 42.048 ±3.263 40.979 ±0.014 41.398 ±4.774 4 20.986 ±0.008 52.922 ±3.524 30.99 ±0.008 53.622 ±4.287 40.988 ±0.009 52.84 ±5.344 Average 0.968 ±0.044 40.487 ±11.078 obtain candidate optimal measurement locations, and how efficient the system is because is the generated paths are all close to the different initial locations of the ASVs. It has been observed that evidently, with a larger number of ASVs, a larger number of measurements will be performed, this will inevitably yield better R2 scores. The obtained GPs manage to correctly fit the data (4 ASVs almost always obtain R2=0.99 ) and approximate the real function with ease. Furthermore, the Voronoi Diagram really enhances robustness and provides scalability. Without the VD, the vehicles will perform measurements on the same regions, and consequently, provide redundant information. With the proposed VD system, the fleet cooperates unintentionally while monitoring. Table 5.16 shows the summarized results of the methods. MOO using GA presents the best result. In this case, it also presents the average greater number of measurements, but it is not so different from the other methods. MOO seems to be more robust regarding initial 5.5 Multi WQP monitoring using Multi-ASVs 73 Table 5.14 Multi Objective Optimization parameter values. Parameter Value Optimizer NSGA-II [86] Individual 2D Gene: (x, y) Population 50 Generations 150 Selection Binary Tournament size =2 Crossover Simulated Binary η=15 ; probcx =0.9 Mutation Polynomial Mutation η=20 ; probmut =0.01 Table 5.15 Results for MOO strategy using GA for different number of ASVs and different number of sensors available. # ASVs # SAvg. R2Avg. M 2 20.953 ±0.027 28.049 ±2.991 30.949 ±0.049 28.774 ±2.894 40.963 ±0.027 29.273 ±3.513 3 20.979 ±0.041 40.045 ±5.002 30.984 ±0.014 42.747 ±4.646 40.987 ±0.01 43.122 ±4.713 4 20.99 ±0.01 53.053 ±8.996 30.989 ±0.014 52.873 ±10.924 40.987 ±0.018 51.722 ±12.033 Average 0.976 ±0.031 41.073 ±12.055 conditions, presenting a lower standard deviation value of R2 . In that sense, the proposed method should include the usage of MOO using GA for accomplishing monitoring of multiple WQPs through multiple ASVs using the Bayesian optimization framework. Table 5.16 Summary of the results using Multi-Water Quality Parameter strategies. Strategy Avg. R2Avg. M MFE Coupled 0.969 ±0.039 40.628 ±10.948 MFE Decoupled 0.968 ±0.044 40.487 ±11.078 MOO using GA 0.976 ±0.031 41.073 ±12.055 74 Chapter 5. Results 0 2000 4000 6000 8000 x [m] 0 2000 4000 6000 8000 10000 12000 14000 y [m] Pareto Ef. Front Other Measurements Other vehicles Current vehicle Figure 5.6 Example of the Multi-ASV Multi-WQP system. A Pareto Efficient set is generated according to the current values returned by the selected acquisition function and the constraint function. The next measurement location is selected to be the closest to the current location of the vehicle and is adapted according to the adaptation function. 5.6 Comparison with other methods Before the results were discussed, a comparison with three different environmental monitoring approaches was performed. The methods are: the Predictive Entropy Search for Multi-Objective Bayesian optimization with Constraints (PESMOC) [ 71 ], applied to the 5.6 Comparison with other methods 75 monitoring objective, the TSP-Genetic Algorithm for monitoring [ 56 ], and finally a lawnmower efficient method for monitoring. The three approaches are modified if possible to achieve maximum exploration and are tested in the same scenario (including ASVs constraints) as the proposed method. Evidently, the proposed system (Figure 3.6) is the same, with the difference in the evaluation component of the acquisition function. The three methods are defined as follows. 1. PESMOC for Environmental Monitoring: In this work, we use the method proposed in [ 71 ] but with some modifications to ensure monitoring. It consists of obtaining the difference between the logarithm of uncertainty of a predictive distribution (PD) σPD i and the average of logarithms of uncertainties of Conditioned PD σCPD i (conditioned to x|X∗(m) ). X∗(m) is one of the m different locations of a supposed Pareto set. For a complete explanation of PESMOC, see [ 71 ]. The PES expression is taken directly from the work, considering deterministic constraints, and has the form of α≈ S ∑ s=1logσPD s(x)−1 MS MS ∑ m=1 logσCPD i(x|X∗(m))(5.2) This expression can be evaluated with the coupled or decoupled strategy. As shown in the expression above, the coupled evaluation sums up the differences for each AF objective αs . The decoupled version considers one difference at a time.The parameters that suffered changes in order to fit the purpose of exploration are as follows. • Pareto set X∗ : Since the objective of this work can be thought of as minimizing uncertainty, the Pareto set is taken as the positions where the sum of the predicted standard deviations reaches its maximum values. • Conditioning p(y|D,x,X∗) : Conditioning is carried out through a cloned GP model to include a supposed evaluation according to the items of the Pareto set. • Monte Carlo Sampling MS : For efficient evaluations, only one point of the Pareto set is used. Therefore, the number of Monte Carlo samples is reduced to one. Indeed, this sacrifices accuracy but definitely improves computational efficiency, which has been observed to be less efficient than our method because it needs to calculate the GP Regression twice for each water quality parameter or objective. Voronoi regions were used to define the constraints, so that each vehicle could only obtain a PESMOC value within its coverage region. With M=1 , as predicted, and using the same simulation sensors as in the evaluation of the proposed method with the best parameter, λ=0.5the results are shown in Table 5.17. 2. TSP-based Environmental Monitoring: In [ 56 ], a set of 60 waypoints were defined on the shore of Ypacarai Lake. Subsequently, the best TSP solution (waypoint visit order) was found by a GA evolved to 76 Chapter 5. Results Table 5.17 Results for PESMOC strategy for different number of ASVs and different number of sensors available. # ASVs # SAvg. R2Avg. M 2 20.907 ±0.056 26.824 ±4.727 30.91 ±0.057 27.174 ±4.358 40.918 ±0.066 28.021 ±4.213 3 20.983 ±0.01 41.82 ±4.161 30.984 ±0.01 43.02 ±4.49 40.983 ±0.008 42.721 ±4.612 4 20.995 ±0.002 54.232 ±5.024 30.995 ±0.002 54.662 ±3.851 40.995 ±0.003 55.469 ±3.168 Average 0.962 ±0.052 41.269 ±12.137 optimize the exploration of Ypacarai Lake. This method was discussed in Chapter 2. Contrary to the continuous measurement approach stated in the aforementioned work [ 56 ], for comparison, the monitoring system is modified so that the vehicle can take measurements only while the ASV is not moving. The distance between the measurement locations is the same as that proposed in this work, which is based on the length scale. Therefore, the ASV travels from point to point, making measurements every l=λ×min{`s} meters. Whenever the total distance traveled reaches 15,000[m], the ASV stops, performs a final measurement, and the mission ends. To account for multiple vehicles, the available shoreline points were divided into groups according to the previously obtained visit order. This is the only difference regarding the usage of multiple ASVs w.r.t. the original proposed system. The results are shown in Table 5.18. 3. Lawnmower applied for monitoring/exploration: The classic lawnmower implementation is also a straightforward approach that seems to be efficient for monitoring, due to its definition. For comparison, the lawnmower strategy found in [ 93 ] was used as part of the acquisition function component, obtaining straight paths that will intend to cover the water body in this problem. The generation of paths considers the Voronoi regions created at execution time and generates waypoint locations that are equally spaced according to the expected prior length scale of the GP. The classic implementation can only deliver paths that are left-to-right and bottom-up (or vice versa, for both pairs of directions), so the vehicles necessarily need to initially be location on the far-most edge of the water body. The results can then only be obtained for multiple sets of simulations that vary the WQP maps but not the initial vehicle locations. The results are available in Table 5.19. None of the comparison methods performed on average better than the proposed method. In Table 5.20, a summary of the comparison methods can be observed, and none of them 5.6 Comparison with other methods 77 Table 5.18 Results for TSP-based GA strategy for different number of ASVs and different number of sensors available. # ASVs # SAvg. R2Avg. M 2 20.837 ±0.086 32.522 ±3.344 30.848 ±0.082 33.574 ±3.065 40.853 ±0.082 34.722 ±2.923 3 20.891 ±0.107 50.171 ±4.017 30.902 ±0.093 51.223 ±3.636 40.903 ±0.093 51.898 ±3.491 4 20.962 ±0.019 60.2±0.4 30.96 ±0.023 60.1±0.944 40.958 ±0.022 60.125 ±0.331 Average 0.902 ±0.088 48.281 ±11.39 Table 5.19 Results for Lawnmower strategy for different number of ASVs and different number of sensors available. # ASVs # SAvg. R2Avg. M 2 20.946 ±0.036 28.9±2.447 30.957 ±0.023 29.75 ±2.808 40.953 ±0.024 30.2±2.421 3 20.815 ±0.109 25.15 ±3.864 30.848 ±0.089 26.25 ±3.819 40.853 ±0.088 26.9±2.488 4 20.739 ±0.132 16.75 ±4.448 30.779 ±0.105 18.25 ±5.29 40.791 ±0.089 18.5±5.045 Average 0.853 ±0.115 24.517 ±6.258 shows better behavior than any of the methods shown in Table 5.16, which contains the results of the proposed method. Of course, because lawnmower and TSP-based using GA are offline methods, they cannot obtain similar R2 values. PESMOC, on the other hand, is an online method that reached a satisfactory average value but is not as good as the proposed method. Figure 5.7 shows a comparison between the methods with respect to the distance of travel. To obtain this figure, the mean values of R2 and the standard deviations were obtained considering the different measurements performed for each method during each simulation. It can be regarded as the average R2 per distance obtained using the different methods. It shows that the MOO using the GA method is not only at the end of the experiments, but is better during the missions because it always presents a better coefficient of determination. Not that the TSP-based method and the lawmower methods have different starting values because their initial positions are selected according to the 78 Chapter 5. Results generated offline paths rather than random initial locations at the shore of the water body. Table 5.20 Summary of the comparison method results. Strategy Avg. R2Avg. M PESMOC 0.962 ±0.052 41.269 ±12.137 TSP-based GA 0.902 ±0.088 48.281 ±11.39 Lawnmower 0.853 ±0.115 24.517 ±6.258 0 2500 5000 7500 10000 12500 15000 Distance [m] 0.2 0.0 0.2 0.4 0.6 0.8 1.0 R 2( x ) Proposed Approach PESMOC TSP-based GA Lawnmower Figure 5.7 Average interpolated R2 value obtained per distance for the proposed approach, as well as the comparison approaches applied to online monitoring. 5.7 Summary of the Results After the obtained results were observed, analyzed, and discussed, a complete summary of them has been done and is as following: • Gaussian Process can be used as surrogate models, since they can efficiently provide regressions using a limited amount of data. With our proposal, as well as with the others, the average coefficient of determination is close to 1, implying that the selected model is very good when approximating unknown models. • Voronoi Regions provide a good systematic approach for multiagent cooperation. This approach only allows for a scalable and robust region partitioning and additionally, in the case of our approach, it can be updated to consider only available ASVs, which is very helpful whenever some agents can be added or removed from the system, in the middle of a mission. • The Multi-Objective problem is appropriately addressed, and is, to best knowledge of the author, one of the firsts works to perform intelligent online water quality 5.7 Summary of the Results 79 parameter measuring for model obtainment. This work expands the Bayesian optimization method to account for multiple objectives (multiple acquisition/utility functions). • Efficiency of Multi-Function Estimation: MFE using multiple vehicles is a quick method for solving multi-objective problems such as this case. Provides the same level of performance as the MOO approach. • Efficiency of Multi-Objective Optimization: It is one of the first approaches that intelligently explore a large-scale environment in the search of data acquisition. The path, which are generated during the mission, is adapted to any scenario and considers both: the acquired data and the endurance of the agent. According to the performed simulations, the MOO approach showed the best results, not only considering the average final R2S obtained, but also when considering the traveled distance. • Comparison with PESMOC: PESMOC is the online method that mos closely resembles the proposed system in this thesis, but the results are not as good. Insights of the method show that it can be trapped in considering points that are too far away, failing to obtain the surrounding behavior but only the behavior between a starting location and the other edges of the map. • Comparison with Genetic Algorithm TSP-based coverage approach: Focusing on the fact that for real applications the measurements should be performed in a discrete manner, the method proposed in this work performs better. As in the LM, this algorithm presents offline paths, so the same differences in efficiency can also be concluded. • Comparison with Lawnmower: Specifically, with this monitoring approach, our method focuses on efficient data acquisition, rather than global, offline, geometricbased data acquisition. The LM method does not account for whether the obtained data fit better or worse to the underlying surrogate model. 86 Bibliography [9] M. Arzamendia, D. Gregor, D. G. Reina, and S. Toral, “An evolutionary approach to constrained path planning of an autonomous surface vehicle for maximizing the covered area of ypacarai lake,” Soft Computing, vol. 23, no. 5, pp. 1723–1734, 2019. [10] S. Y. Luis, D. G. Reina, and S. T. Marín, “A deep reinforcement learning approach for the patrolling problem of water resources through autonomous surface vehicles: The ypacarai lake case,” IEEE Access, vol. 8, pp. 204 076–204 093, 2020. [11] S. Y. Luis, D. G. Reina, and S. L. Toral, “A multiagent deep reinforcement learning approach for path planning in autonomous surface vehicles: The ypacaraí lake patrolling case,” IEEE Access, vol. 9, pp. 17 084–17 099, 2021. [12] C. Chauvin-Hameau, “Informative path planning for algae farm surveying,” 2020. [13] Z. Peng, J. Wang, D. Wang, and Q.-L. Han, “An overview of recent advances in coordinated control of multiple autonomous surface vehicles,” IEEE Transactions on Industrial Informatics, vol. 17, no. 2, pp. 732–745, 2020. [14] M. Arzamendia, D. Gregor, D. G. Reina, and S. Toral, “A path planning approach of an autonomous surface vehicle for water quality monitoring using evolutionary computation,” in Technology for Smart Futures. Springer, 2018, pp. 55–73. [15] M. Arzamendia, I. Espartza, D. G. Reina, S. Toral, and D. Gregor, “Comparison of eulerian and hamiltonian circuits for evolutionary-based path planning of an autonomous surface vehicle for monitoring ypacarai lake,” Journal of Ambient Intelligence and Humanized Computing, vol. 10, no. 4, pp. 1495–1507, 2019. [16] N. Karapetyan, A. Braude, J. Moulton, J. A. Burstein, S. White, J. M. O’Kane, and I. Rekleitis, “Riverine coverage with an autonomous surface vehicle over known environments,” in 2019 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), 2019, pp. 3098–3104. [17] G. Hitz, E. Galceran, M. Garneau, F. Pomerleau, and R. Siegwart, “Adaptive continuous-space informative path planning for online environmental monitoring,” Journal of Field Robotics, vol. 34, no. 8, pp. 1427–1449, 2017. [18] L. Jin, J. Rückin, S. H. Kiss, T. Vidal-Calleja, and M. Popović, “Adaptive-resolution gaussian process mapping for efficient uav-based terrain monitoring,” arXiv preprint arXiv:2109.14257, 2021. [19] P. Morere, R. Marchant, and F. Ramos, “Sequential bayesian optimization as a pomdp for environment monitoring with uavs,” in 2017 IEEE International Conference on Robotics and Automation (ICRA). IEEE, 2017, pp. 6381–6388. [20] C. Williams and C. Rasmussen, Gaussian processes for machine learning. MIT press Cambridge, MA, 2006, vol. 2. Bibliography 87 [21] S. Kemna, J. G. Rogers, C. Nieto-Granda, S. Young, and G. S. Sukhatme, “Multirobot coordination through dynamic voronoi partitioning for informative adaptive sampling in communication-constrained environments,” in 2017 IEEE International Conference on Robotics and Automation (ICRA). IEEE, 2017, pp. 2124–2130. [22] Z. Liu, Y. Zhang, X. Yu, and C. Yuan, “Unmanned surface vehicles: An overview of developments and challenges,” Annual Reviews in Control, vol. 41, pp. 71–93, 2016. [23] C. Zhou, S. Gu, Y. Wen, Z. Du, C. Xiao, L. Huang, and M. Zhu, “The review unmanned surface vehicle path planning: Based on multi-modality constraint,” Ocean Engineering, vol. 200, p. 107043, 2020. [24] L. Bottarelli, M. Bicego, and A. F. J. Blum, “Orienteering-based informative path planning for environmental monitoring,” Engineering Applications of Artificial Intelligence, vol. 77, pp. 46–58, 2019. [25] L. Iftekhar, H. N. Rahman, and I. Rahman, “Area coverage algorithms for networked multi-robot systems,” in New Developments and Advances in Robot Control. Springer, 2019, pp. 301–320. [26] Z. Zhen, D. Xing, and C. Gao, “Cooperative search-attack mission planning for multi-uav based on intelligent self-organized algorithm,” Aerospace Science and Technology, vol. 76, pp. 402–411, 2018. [27] C. Sampedro, H. Bavle, J. L. Sanchez-Lopez, R. A. S. Fernández, A. Rodríguez- Ramos, M. Molina, and P. Campoy, “A flexible and dynamic mission planning architecture for uav swarm coordination,” in 2016 International Conference on Unmanned Aircraft Systems (ICUAS). IEEE, 2016, pp. 355–363. [28] T. T. Mac, C. Copot, D. T. Tran, and R. De Keyser, “Heuristic approaches in robot path planning: A survey,” Robotics and Autonomous Systems, vol. 86, pp. 13–28, 2016. [29] H. Shin and J. Chae, “A performance review of collision-free path planning algorithms,” Electronics, vol. 9, no. 2, p. 316, 2020. [30] I. Chaari, A. Koubaa, H. Bennaceur, A. Ammar, M. Alajlan, and H. Youssef, “Design and performance analysis of global path planning techniques for autonomous mobile robots in grid environments,” International Journal of Advanced Robotic Systems, vol. 14, no. 2, p. 1729881416663663, 2017. [31] C. Zammit and E.-J. Van Kampen, “Comparison between a* and rrt algorithms for uav path planning,” in 2018 AIAA guidance, navigation, and control conference, 2018, p. 1846. [32] B. Paden, M. Čáp, S. Z. Yong, D. Yershov, and E. Frazzoli, “A survey of motion planning and control techniques for self-driving urban vehicles,” IEEE Transactions on intelligent vehicles, vol. 1, no. 1, pp. 33–55, 2016. 88 Bibliography [33] H. Abraham, B. Seppelt, B. Mehler, and B. Reimer, “What’s in a name: Vehicle technology branding & consumer expectations for automation,” in Proceedings of the 9th international conference on automotive user interfaces and interactive vehicular applications, 2017, pp. 226–234. [34] C. Wickens, “Automation stages & levels, 20 years after,” Journal of Cognitive Engineering and Decision Making, vol. 12, no. 1, pp. 35–41, 2018. [35] P. Tokekar, J. Vander Hook, D. Mulla, and V. Isler, “Sensor planning for a symbiotic uav and ugv system for precision agriculture,” IEEE Transactions on Robotics, vol. 32, no. 6, pp. 1498–1511, 2016. [36] C. Liu, Q. Mao, X. Chu, and S. Xie, “An improved a-star algorithm considering water current, traffic separation and berthing for vessel path planning,” Applied Sciences, vol. 9, no. 6, p. 1057, 2019. [37] L. Bayındır, “A review of swarm robotics tasks,” Neurocomputing, vol. 172, pp. 292–321, 2016. [38] A. Kapoutsis, “Towards a fully autonomous and cooperative deployment of multirobot teams for exploration and coverage in unknown or partially known environments,” Ph.D. dissertation, ∆ηµ o κρ´ ιτει o Πανεπιστ ´ ηµι o Θρ´ ακης ( ∆ΠΘ ). Σχoλ´ ηΠoλυτεχνικ ´ η. Tµ´ ηµα Hλεκτρoλóγων ..., 2017. [39] D. Brandtner and M. Saska, “Coherent swarming of unmanned micro aerial vehicles with minimum computational and communication requirements,” in 2017 European Conference on Mobile Robots (ECMR). IEEE, 2017, pp. 1–6. [40] J. Alonso-Mora, S. Baker, and D. Rus, “Multi-robot formation control and object transport in dynamic environments via constrained optimization,” The International Journal of Robotics Research, vol. 36, no. 9, pp. 1000–1021, 2017. [41] T. Yang, S. Hsiung, K. P. CH. Kuo, YD. Tsai, Y. Hsieh, Z. Shen, J. Feng, and C. Kuo, “Development of unmanned surface vehicle for water quality monitoring and measurement,” in 2018 IEEE International Conference on Applied System Invention (ICASI). IEEE, 2018, pp. 566–569. [42] J. Paez, J. Villa, J. Cabrera, and E. Yime, “Implementation of an unmanned surface vehicle for environmental monitoring applications,” in 2018 IEEE 2nd Colombian Conference on Robotics and Automation (CCRA). IEEE, 2018, pp. 1–6. [43] J. Sánchez-García, J. García-Campos, M. Arzamendia, D. G. Reina, S. Toral, and D. Gregor, “A survey on unmanned aerial and aquatic vehicle multi-hop networks: Wireless communications, evaluation tools and applications,” Computer Communications, vol. 119, pp. 43–65, 2018. [44] V. A. Jorge, R. Granada, R. G. Maidana, D. A. Jurak, G. Heck, A. P. Negreiros, D. H. Dos Santos, L. M. Gonçalves, and A. M. Amory, “A survey on unmanned surface vehicles for disaster robotics: Main challenges and directions,” Sensors, vol. 19, no. 3, p. 702, 2019. Bibliography 89 [45] Q. Yang, S.-J. Jang, and S.-J. Yoo, “Q-learning-based fuzzy logic for multi-objective routing algorithm in flying ad hoc networks,” Wireless Personal Communications, pp. 1–24, 2020. [46] D. Angley, B. Ristic, W. Moran, and B. Himed, “Search for targets in a risky environment using multi-objective optimisation,” IET Radar, Sonar & Navigation, vol. 13, no. 1, pp. 123–127, 2018. [47] S. Sawadsitang, D. Niyato, P. S. Tan, P. Wang, and S. Nutanong, “Multi-objective optimization for drone delivery,” in 2019 IEEE 90th Vehicular Technology Conference (VTC2019-Fall). IEEE, 2019, pp. 1–5. [48] K. I. Kilic, O. Gemikonakli, and L. Mostarda, “Multi-objective priority based heuristic optimization for region coverage with uavs,” in International Conference on Advanced Information Networking and Applications. Springer, 2020, pp. 768–779. [49] Z. Gosiewski and K. Kwaśniewski, “Time minimization of rescue action realized by an autonomous vehicle,” Electronics, vol. 9, no. 12, p. 2099, 2020. [50] C. Ju and H. I. Son, “Multiple uav systems for agricultural applications: control, implementation, and evaluation,” Electronics, vol. 7, no. 9, p. 162, 2018. [51] R. Almadhoun, T. Taha, L. Seneviratne, and Y. Zweiri, “A survey on multi-robot coverage path planning for model reconstruction and mapping,” SN Applied Sciences, vol. 1, no. 8, p. 847, 2019. [52] K. Yu, J. O’Kane, and P. Tokekar, “Coverage of an environment using energyconstrained unmanned aerial vehicles,” in 2019 International Conference on Robotics and Automation (ICRA). IEEE, 2019, pp. 3259–3265. [53] L. Schmid, M. Pantic, R. Khanna, L. Ott, R. Siegwart, and J. Nieto, “An efficient sampling-based method for online informative path planning in unknown environments,” IEEE Robotics and Automation Letters, vol. 5, no. 2, pp. 1500–1507, 2020. [54] F. Peralta, M. Arzamendia, D. Gregor, D. G. Reina, and S. Toral, “A comparison of local path planning techniques of autonomous surface vehicles for monitoring applications: The ypacarai lake case-study,” Sensors, vol. 20, no. 5, p. 1488, 2020. [55] Itaipú, “Centro internacional de hidroinformática,” https://hidroinformatica.itaipu. gov.py/, 2016, accessed: 2020-10-15. [56] M. Arzamendia, D. Gutierrez, S. Toral, D. Gregor, E. Asimakopoulou, and N. Bessis, “Intelligent online learning strategy for an autonomous surface vehicle in lake environments using evolutionary computation,” IEEE Intelligent Transportation Systems Magazine, vol. 11, no. 4, pp. 110–125, 2019. 90 Bibliography [57] S. Yanes, D. Gutiérrez-Reina, and S. Toral, “A dimensional comparison between evolutionary algorithm and deep reinforcement learning methodologies for autonomous surface vehicles with water quality sensors,” Sensors, vol. 21, no. 8, p. 2862, 2021. [58] M. J. T. Kathen, I. J. Flores, and D. G. Reina, “An informative path planner for a swarm of asvs based on an enhanced pso with gaussian surrogate model components intended for water monitoring applications,” Electronics, vol. 10, no. 13, p. 1605, 2021. [59] M. J. Ten Kathen, I. J. Flores, D. G. Reina, and A. T. Córdoba, “Autonomous monitoring system for water resources based on pso and gaussian process,” in 2021 IEEE Congress on Evolutionary Computation (CEC). IEEE, 2021, pp. 1777–1784. [60] M. J. Ten Kathen, I. J. Flores, and D. G. Reina, “A comparison of pso-based informative path planners for autonomous surface vehicles for water resource monitoring,” in 7th International Conference on Machine Learning Technologies (ICMLT) (ICMLT 2022). ACM, 2022 in press. [61] M. J. Ten Kathen, D. G. Reina, and I. J. Flores, “A comparison of pso-based informative path planners for detecting pollution peaks of the ypacarai lake with autonomous surface vehicles,” in International Conference on Optimization and Learning OLA’2022, 2022 in press. [62] S. Y. Luis, F. Peralta, A. T. Córdoba, Á. R. del Nozal, S. T. Marín, and D. G. Reina, “An evolutionary multi-objective path planning of a fleet of asvs for patrolling water resources,” Engineering Applications of Artificial Intelligence, vol. 112, p. 104852, 2022. [63] M. Popović, T. Vidal-Calleja, G. Hitz, J. Chung, I. Sa, R. Siegwart, and J. Nieto, “An informative path planning framework for uav-based terrain monitoring,” Autonomous Robots, pp. 1–23, 2020. [64] Y. Qiming, Z. Jiandong, and S. Guoqing, “Modeling of uav path planning based on imm under pomdp framework,” Journal of Systems Engineering and Electronics, vol. 30, no. 3, pp. 545–554, 2019. [65] R. Marchant and F. Ramos, “Bayesian optimisation for informative continuous path planning,” in 2014 IEEE International Conference on Robotics and Automation (ICRA). IEEE, 2014, pp. 6136–6143. [66] R. Martinez-Cantin, N. D. Freitas, E. Brochu, J. Castellanos, and A. Doucet, “A bayesian exploration-exploitation approach for optimal online sensing and planning with a visually guided mobile robot,” Autonomous Robots, vol. 27, no. 2, pp. 93–103, 2009. [67] R. Oliveira, L. Ott, V. Guizilini, and F. Ramos, “Bayesian optimisation for safe navigation under localisation uncertainty,” in Robotics Research. Springer, 2020, pp. 489–504. Bibliography 91 [68] J. Gao, L. Zeng, C. Cao, W. Ye, and X. Zhang, “Multi-objective optimization for sensor placement against suddenly released contaminant in air duct system,” in Building simulation, vol. 11. Springer, 2018, pp. 139–153. [69] S. Pourshahabi, M. R. Nikoo, E. Raei, and J. F. Adamowski, “An entropy-based approach to fuzzy multi-objective optimization of reservoir water quality monitoring networks considering uncertainties,” Water Resources Management, vol. 32, no. 13, pp. 4425–4443, 2018. [70] D. Reina, H. Tawfik, and S. Toral, “Multi-subpopulation evolutionary algorithms for coverage deployment of uav-networks,” Ad Hoc Networks, vol. 68, pp. 16–32, 2018. [71] E. C. Garrido-Merchán and D. Hernández-Lobato, “Predictive entropy search for multi-objective bayesian optimization with constraints,” Neurocomputing, vol. 361, pp. 50–68, 2019. [72] H. M. Torun, M. Swaminathan, A. K. Davis, and M. L. F. Bellaredj, “A global bayesian optimization algorithm and its application to integrated system design,” IEEE Transactions on Very Large Scale Integration (VLSI) Systems, vol. 26, no. 4, pp. 792–802, 2018. [73] N. Karapetyan, J. Moulton, J. Lewis, A. Li, A. Quattrini, J. O’Kane, and I. Rekleitis, “Multi-robot dubins coverage with autonomous surface vehicles,” in 2018 IEEE International Conference on Robotics and Automation (ICRA). IEEE, 2018, pp. 2373–2379. [74] A. Vasilijević, ÐNađ., F. Mandić, N. Mišković, and Z. Vukić, “Coordinated navigation of surface and underwater marine robotic vehicles for ocean sampling and environmental monitoring,” IEEE/ASME transactions on mechatronics, vol. 22, no. 3, pp. 1174–1184, 2017. [75] A. Dutta, A. Ghosh, and O. P. Kreidl, “Multi-robot informative path planning with continuous connectivity constraints,” in 2019 International Conference on Robotics and Automation (ICRA), 2019, pp. 3245–3251. [76] E. Camci, D. Kripalani, L. Ma, E. Kayacan, and M. Khanesar, “An aerial robot for rice farm quality inspection with type-2 fuzzy neural networks tuned by particle swarm optimization-sliding mode control hybrid algorithm,” Swarm and evolutionary computation, vol. 41, pp. 1–8, 2018. [77] T. Gao, H. Emadi, H. Saha, J. Zhang, A. Lofquist, A. Singh, B. Ganapathysubramanian, S. Sarkar, A. K. Singh, and S. Bhattacharya, “A novel multirobot system for plant phenotyping,” Robotics, vol. 7, no. 4, p. 61, 2018. [78] D. P. Nicholson, A. P. M. Michel, S. D. Wankel, K. Manganini, R. A. Sugrue, Z. O. Sandwith, and S. A. Monk, “Rapid mapping of dissolved methane and carbon dioxide in coastal ecosystems using the chemyak autonomous surface vehicle,” Environmental science & technology, vol. 52, no. 22, pp. 13 314–13 324, 2018. 92 Bibliography [79] F. Archetti and A. Candelieri, Bayesian Optimization and Data Science. Springer, 2019. [80] T. Abbasi and S. Abbasi, “Conventional’indices for determining fitness of waters for different uses,” Water Quality Indices; Elsevier: Amsterdam, The Netherlands, pp. 25–62, 2012. [81] M. Buisson-Fenet, F. Solowjow, and S. Trimpe, “Actively learning gaussian process dynamics,” in Learning for dynamics and control. PMLR, 2020, pp. 5–15. [82] J. Görtler, R. Kehlbeck, and O. Deussen, “A visual exploration of gaussian processes,” Distill, 2019, https://distill.pub/2019/visual-exploration-gaussian- processes. [83] C. F. Falk and M. Muthukrishna, “Parsimony in model selection: Tools for assessing fit propensity.” Psychological Methods, 2021. [84] D. C. Liu and J. Nocedal, “On the limited memory bfgs method for large scale optimization,” Mathematical programming, vol. 45, no. 1, pp. 503–528, 1989. [85] U. Noè and D. Husmeier, “On a new improvement-based acquisition function for bayesian optimization,” arXiv preprint arXiv:1808.06918, 2018. [86] K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan, “A fast and elitist multiobjective genetic algorithm: Nsga-ii,” IEEE transactions on evolutionary computation, vol. 6, no. 2, pp. 182–197, 2002. [87] J. Blank and K. Deb, “Pymoo: Multi-objective optimization in python,” IEEE Access, vol. 8, pp. 89 497–89509, 2020. [88] Emlid, “Emlid reach docs,” https://docs.emlid.com/reach/, 2019, [Online; accessed 17-05-2022]. [89] R. P. Foundation, “Raspberyy pi 4 model b,” https://www.raspberrypi.com/products/ raspberry-pi-4-model-b/, 2019, [Online; accessed 17-05-2022]. [90] Nvidia, “Nvidia jetson agx xavier,” nvidia.com/es-es/autonomous-machines/ embedded-systems/jetson-agx-xavier/, 2019, [Online; accessed 17-05-2022]. [91] T. A. Morel, J. M. Manzano, G. Bejarano, L. Orihuela et al., “Modelling and identification of an autonomous surface vehicle: Technical report,” Informes Técnicos, 2022. [92] Libelium, “Smart water sensor guide,” https://development.libelium.com/smart_ water_sensor_guide/sensors, 2019, [Online; accessed 17-05-2022]. [93] S. A, D. Ingram, J. Dinius, K. Chawla, A. Raffin, and A. Paques, “Pythonrobotics: a python code collection of robotics algorithms,” arXiv preprint arXiv:1808.10703, 2018. Bibliography 93 [94] F. Peralta, M. Arzamendia, D. Gregor, K. Cikel, M. Santacruz, D. G. Reina, and S. Toral, “Development of a simulator for the study of path planning of an autonomous surface vehicle in lake environments,” in 2019 IEEE CHILEAN Conference on Electrical, Electronics Engineering, Information and Communication Technologies (CHILECON). IEEE, 2019, pp. 1–6. [95] F. P. Samaniego, D. G. Reina, S. L. T. Marín, M. Arzamendia, and D. O. Gregor, “A bayesian optimization approach for water resources monitoring through an autonomous surface vehicle: The ypacarai lake case study,” IEEE Access, vol. 9, pp. 9163–9179, 2021. [96] F. Peralta, D. G. Reina, S. Toral, M. Arzamendia, and D. Gregor, “A bayesian optimization approach for multi-function estimation for environmental monitoring using an autonomous surface vehicle: Ypacarai lake case study,” Electronics, vol. 10, no. 8, p. 963, 2021. [97] F. Peralta, D. G. Reina, and S. Toral, “Water quality online modeling using multiobjective and multi-agent bayesian optimization with region partitioning,” in Engineering Applications of Artificial Intelligence. Elsevier, 2022 under review. [98] F. Peralta, S. Yanes, D. G. Reina, and S. Toral, “Monitoring water resources through a bayesian optimization-based approach using multiple surface vehicles: The ypacarai lake case study,” in 2021 IEEE Congress on Evolutionary Computation (CEC). IEEE, 2021, pp. 1511–1518. [99] F. Peralta, D. G. Reina, and S. Toral, “Towards an online water quality monitoring system of dynamic environments using an autonomous surface vehicle,” in International Conference on Optimization and Learning OLA’2022, 2022 in press. [100] F. Peralta, M. Pearce, M. Poloczek, D. G. Reina, S. Toral, and J. Branke, “Multiobjective path planning for environmental monitoring using an autonomous surface vehicle,” Genetic and Evolutionary Computation Conference Companion (GECCO ’22 Companion), 2022.