Multi-robot sensor networks for thermosolar power plants

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Tesis Doctoral Ingeniería Industrial Multi-Robot Sensor Networks for Thermosolar Power Plants Autor: Javier García Martín Directores: Jose Maria Maestre Torreblanca Eduardo Fernandez Camacho Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2023 Tesis Doctoral Ingeniería Industrial Multi-Robot Sensor Networks for Thermosolar Power Plants Autor: Javier García Martín Directores: Jose Maria Maestre Torreblanca Catedrático de Universidad Eduardo Fernandez Camacho Catedrático de Universidad Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla 2023 Tesis Doctoral: Multi-Robot Sensor Networks for Thermosolar Power Plants Autor: Javier García Martín Directores: Jose Maria Maestre Torreblanca Eduardo Fernandez Camacho 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 mi familia, mis amigos, y mis maestros Agradecimientos En primer lugar, quiero expresar mi agradecimiento a mis profesores, en especial a Pepe, que hace cuatro años me llamó y me abrió este camino que siempre me había atraído, y a Eduardo, por ser un excelente guía tanto para mí como para el resto de doctorandos. También quiero dar las gracias al resto del grupo de investigación, que ha hecho que trabajar aquí durante estos años haya sido un verdadero placer al haber creado un magnífico ambiente de compañerismo y confianza. También quiero agradecer a los profesores de instituciones extranjeras que me han acogido durante este tiempo, el Dr. Eric Duviella del IMT Nord Europe y el Dr. Takeshi Hatanaka del Tokyo Institute of Technology, por la hospitalidad que me han brindado. Como no puede ser de otro modo quiero expresar también mi más sincero agradecimiento a mis padres por su constante apoyo y por haberme brindado todo lo que necesitaba durante mis años de doctorado. Gracias por su incondicional amor y por haberme dado la oportunidad de recibir una excelente educación desde que nací. Sin ellos, no habría podido lograr todo lo que he conseguido hasta ahora. Su compromiso, su esfuerzo y su dedicación son algo que nunca podré pagar. Gracias por ser mis pilares más fuertes y por estar siempre a mi lado, brindándome vuestro apoyo incondicional. Nunca podré agradeceros lo suficiente todo lo que habeis hecho por mí. Del mismo modo, quisiera también rendir homenaje a mis abuelos y abuelas, que siempre trabajaron arduamente para asegurar un futuro mejor para sus descendientes. A pesar de que ellos mismos venían de familias humildes y recibieron poca educación, lucharon para salir adelante, y siempre se esforzaron para que tuvieramos una buena educación sin imaginar nunca que alguno de nosotros pudiera convertirse algún día en doctor. Vuestros nietos siempre estaremos en deuda con vosotros. También quiero expresar mi agradecimiento mas profundo a mi compañera de vida, Inma, por haberme acompañado y apoyado durante esta fase final del doctorado siempre que lo he necesitado. Gracias por tu constante apoyo y amor. Eres mi hogar. No puedo olvidarme de mi familia de Urbis, Cesar, Jesús y Santi, por todo el acompañamiento y el apoyo que me han brindado durante estos años. Vivir con vosotros ha sido una experiencia II Contents 3.2.1 Mixed Integer Linear Programming Problem Formulation 29 3.2.2 Tuning Costs 30 3.3 Simulation Environment 31 3.4 Results 33 3.5 Conclusions 39 4 BaB and GA for TA MRTA Problem in the PTC-CSP 41 4.1 Variables 42 4.2 Cost Function 43 4.3 Branch and Bound Algorithm 47 4.4 Genetic Algorithm 51 4.5 Initial Solutions 54 4.6 Case Study 56 4.6.1 Scenarios 56 4.6.2 Genetic Algorithm tuning 56 4.7 Results 57 4.7.1 Scenario I 57 4.7.2 Scenario II 59 4.7.3 Scenario III 59 4.7.4 Scenario IV 59 4.7.5 Performance Comparison 60 4.7.6 Computational Costs 61 4.7.7 Montecarlo Studio 62 4.8 Conclusions 63 5 Game Theory Based Clustering for MRTA 65 5.1 Problem Statement 66 5.1.1 MRTA Objective Function 67 5.1.2 Cluster Evaluation 67 5.1.3 Clustering Evaluation 69 5.1.4 Proposed Control Scheme 69 5.2 Game Theory Viewpoint 71 5.2.1 Estimating the Shapley Value 73 5.3 Clustering Algorithm 76 5.4 Case Study 79 5.4.1 Centralized Clustering 79 5.4.2 SVCA Clustering 79 5.4.3 Genetic Algorithms (GA)-Based Clustering 88 5.4.4 Averaging Results 94 5.4.5 Scalability Analisys 94 5.5 Conclusions 95 6 Predictive MRTA with dynamic tasks 97 Contents III 6.1 Problem Statement 97 6.1.1 Problem Formulation 104 6.2 Predictive Multi-Robot Task Allocation Algorithm (PMRTA) 105 6.2.1 Predictive Receding-Horizon Multi-Robot Task Allocation Algorithm 107 6.3 Simulations and Results 108 6.3.1 High-Fidelity Simulation 111 6.4 Conclusions 117 Part III Direct Normal Irradiance Estimation 119 7 Bayesian Based Inference Estimation 121 7.1 Problem Statement 121 7.1.1 Information Sources and Processing 123 7.1.2 Outcomes 124 7.2 Task Generation 125 7.3 Case-Study 127 7.3.1 Bayesian estimation parameters 128 7.3.2 Wind and clouds 128 7.4 Results 129 7.5 Conclusions 132 8 Data-Fusion Spatio-Temporal kriging Estimation 135 8.1 Spatial and Temporal DNI mapping using Kriging 136 8.1.1 Kriging 136 8.1.2 Identification of the STV parameters 138 8.1.3 The Point of Min. Semivariance and the Selection of Measurements 139 8.2 Case Study and Results 139 8.2.1 Obtention of the Variogram and Adjustment to Theoretical Model 140 8.2.2 Baseline Variogram and Adjustment 142 8.3 Map Estimation Results 143 8.3.1 Sensor meshes 143 8.3.2 Random instant analysis 144 8.3.3 Averaging results 146 8.4 Conclusions 147 Part IV Conclusions and Future Research 149 9 Conclusions and Future Research 151 9.1 Conclusions 151 9.2 Future Research 153 9.2.1 Robot Fleet Management 153 9.2.2 DNI Estimation and Forecasting 154 IV Contents List of Figures 155 List of Tables 159 Bibliography 161 Glossary ACO Ant Colony Algorithm. 8 ADMM Alternating Direction Method of Multipliers. 153 BaB Branch and Bound. 9, 41, 42, 44, 46–52, 54, 56–63, 156, 159, II, XIII BI Bayesian Inference. 121–123, 152 CD Complex Dependency between tasks. 7 CF Cloud Factor. 20, 121 CSP Concentrated Solar Plant. 3 DCOP Distributed Constraint Optimization Problem. 8 DNI Direct Normal Irradiance. 3, 4, 9, 19–21, 41, 121, 122, 124, 125, 128, 132, 135, 136, 147, 151–153, XII GA Genetic Algorithm. 8, 9, 41, 42, 44, 46, 48, 50–52, 54, 56–63, 66, 88–95, 108–110, 113, 125, 129, 153, 156, 157, 159, 160, II, XV GIS Geographic Information System. 135 HTF Heat Transfer Fluid. 3, 4, 15 IA Instantaneus Assignment. 7, 9, 25, 26, 28, 30, 32, 34, 36, 38, 41, I ID In-schedule Dependency between tasks. 7 LCOE Levelized Costs of Electricity. 3 LP Linear Programming. 97, 105, 152 V VI Glossary MILP Mixed Integer Linear Programming. 8, 25, 26, 28–30, 32, 34, 36, 38, I MPC Model Predictive Control. 4, 97 MR Multi-Robot tasks. 7 MRS Multi-Robot System. 5, 6, 8, 66, 133 MRTA Multi-Robot Task Allocation. 7–9, 25, 26, 28, 30, 32, 34, 36, 38, 39, 41, 42, 44–46, 48, 50, 52, 54, 56, 58, 60, 62, 65–70, 72, 74, 76, 78, 80, 82, 84, 86, 88, 90, 92–97, 105, 117, 125, 129, 132, 152, 159, I, II MRTD Multi-Robot Task Decomposition. 7 MRTP Multi-Robot Task Planning. 6–8 MT Multi-Task robots. 7 ND No Dependency between tasks. 7 NN Nearest neighbor. 54 OAP Optimal Assignment Problem. 8, 25, 30, 33–39, 155 OCONTSOLAR Optimal Control of Thermal Solar Energy Systems. 4, 5, 14, 79, 151, 155 PMRTA Predictive Receding-Horizon Multi-Robot Task Allocation. 107, 109–111, 113, 152, 157, XV PSO Particle Swarm Optimization. 8 PTC-CSP Parabolic Trough Collector Concentrated Solar Plant. 3–5, 9, 15–17, 19, 21, 25, 26, 28, 30, 32, 34, 36, 38, 41, 42, 44, 46, 48, 50, 52, 54, 56, 58, 60, 62, 63, 121, 151, 155, I, II PVP Photovoltaic plant. 3, 136 ROS Robotic Operative System. 111–115, 154, 157 RSN Robotic Sensor Network. 4–6, 25, 121, 151, 154 SOC State of Charge. 14, 38, 39, 104, 153, 155 SR Single-Robot tasks. 7, 9, 66 ST Single-Task robots. 7, 9, 66 STK Spatio-Temporal kriging. 136, 152, 154 STV Spatio-Temporal Variogram. 135–140, 142, 147, 152 SVCA Shapley Value Clustering Algorithm. 76, 77, 79, 81–89, 93–95, 156, 159, 160, II Glossary VII TA Time-Extended Assignment. 7, 9, 41, 42, 44, 46, 48, 50, 52, 54, 56, 58, 60, 62, 66, 97, II TES Thermic Energy Storage. 3, 15 TSP Traveling Salesman Problem. 54 UAV Unmanned Aerial Vehicle. 6, 9, 13–15, 26, 27, 30–35, 38, 41, 47, 56, 59, 63, 112, 121, 122, 127, 151, 155, 159 UGV Unmanned Ground Vehicle. 6, 9, 13–15, 17, 25–27, 30, 31, 33–35, 38, 41, 47, 56, 58, 59, 63, 112, 127, 151, 155, 159 VTOL Vertical Take-Off Landing. 14 WSN Wireless Sensor Network. 4, 154 XD Cross Dependency between tasks. 7, 9 Notation Math RReal numbers CComplex numbers ∥v∥Norm-2 of vector v |A|Determinant of A A⊤Transpose of matrix A A−1Inverse of matrix A InIdentity matrix of size n 1nSquare matrix of ones of size n 1m×nMatrix of ones of size m×n 0nSquare matrix of zeros of size n 0m×nMatrix of zeros of size m×n diag(v)Diagonal matrix from vector v ∂y ∂xPartial derivative of yrespect x Pr(A)Likehood of event A E[X]Expected value of random variable X σ2 XVariance of random variable X ∼fX(x)Distributed following the probability density function fX(x) NµX,σ2 X Gaussian distribution of random variable X , with mean µX and variance σ2 X |A|Cardinal, number of elements in set A OOrder General Environment IX XNotation QEnvironment CSet of cells in which the environment Qis discretized cjCell, element of C qjPosition of the centre of cj CNumber of cells in which Qis discretized NrNumber of rows in which Qis discretized NcNumber of columns in which Qis discretized Robots Set RRobots set RNumber of tasks (cardinal of R) riRobot i, element of set Rwith ri,i=1,2,...,R piPosition of robot ri RGUGVs set giUGV i, element of set RGwith gi,i=1,2,...,RG RGNumber of UGVs (cardinal of RG) RAUAVs set RANumber of UAVs (cardinal of RA) aiUAV i, element of set RAwith ai,i=1,2,...,RA UGV obstacles OGSet of obstacles for UGVs OGNumber of obstacles for UGVs OG jj-esim obstacle for UGVs rG jl l-th vertex of OG j ρG jl Position of the l-th vertex of OG j VGSet of vertixes of UGVs obstacles UAV obstacles OASet of obstacles for UAVs OANumber of obstacles for UAVs OA ii-esim obstacle for UAVs rA il l-th vertex of OA i ρA il Position of the l-th vertex of OA i VASet of vertixes of UAVs obstacles Distance DiWeighted undirected graph associated to ri Notation XI NG iNodes set of graph Diif ri∈RG NA iNodes set of graph Diif ri∈RA EG iEdges set of graph Diif ri∈RG EA iEdges set of graph Diif ri∈RA WG iWeights set of graph Diif ri∈RG WA iWeights set of graph Diif ri∈RA wG i j Element of set WG i wA i j Element of set WA i dDistance function Tasks set TTasks set TNumber of tasks (cardinal of T) tjTask j, element of set Twith tj,j=1,2,...,T Allocation horizon KSet of instants in the allocation horizon kInstant. Element of K KAllocation horizon (cardinal of K) Robot Parameters λiPenalty for moving a certain robot ri viMaximum velocity of robot ri(m/s) vZ iMaximum vertical (land/take off) velocity of robot ri(m/s) wi%of battery consumed by robot riper second of use hi%of battery charge by robot riper second charging biBattery state of charge (%) of robot ri bo iInitial battery state of charge (%) of robot ri bS iMinimum battery state of charge (%) of robot ri Task Parameters ϕjWeight of task tjcompletion time ψjTime required (s) to perform task tjonce it has been started Clouds CF Cloud factor received in a cell (from 0 if there is no cloud to 1if no irradiance is received) Part I Introduction and Problem Setting 1 1 Introduction There is an endless number of things to discover about robotics. A lot of it is just too fantastic for people to believe. Daniel H. Wilson, 2011 Nowadays , climate change is considered one of the top challenges for the future of humanity. For this reason, a significant amount of research has focused on the use of renewable energies during the last decades, decreasing their levelized cost of electricity (LCOE) and reaching lower costs than those of fossil fuel technologies (Timilsina, 2021). Therefore, it is not surprising that the use of solar energy around the world has rapidly increased during the last years. The most common plants for electric production using solar energy are photovoltaic plants (PVPs) and concentrated solar thermal plants (CSPs). While the former have lower LCOE, CSPs can incorporate thermal energy storage (TES) (for example, using molten salt) to generate electricity even during the night, which can be very convenient. CSPs are usually located on large extensions of land with a high solar incidence, being those with parabolic trough collector (PTC-CSP) the most profitable and common type today (Lilliestam and Pitz-Paal, 2018). These plants use parabolic mirrors or solar collectors to concentrate solar energy on pipes located at the focal point of the parabola, where a heat transfer fluid (HTF) circulates absorbing the solar radiation. The hot HTF from the solar collectors is then sent to a power generation plant through a manifold and is used to generate steam to feed a turbine, generating electricity through the Rankine cycle (Camacho and Berenguel, 2012; Camacho et al., 1997). After that, another manifold distributes the cold HTF back to the collectors. In most commercial plants, the solar field is composed of a number of loops connected in parallel and each loop is formed by a number of collectors, usually four, connected in series. Each collector tracks the sun on one axis to maximize the collected energy. PTC-CSPs are usually operated by controlling the total flow that goes through the manifolds, generally with the aim of maintaining the HTF temperature close to a given setpoint. In general, the more DNI received, the greater the flow will have to be so that the HTF temperature 3 4Chapter 1. Introduction does not exceed a certain upper bound and; in the same way, when the DNI is low, the flow will be decreased so that the temperature respects its lower bound. However, this strategy has an important drawback whenever the irradiance is not uniform over the plant. For example, a localized cloud might trigger a flow decrease, rising HTF temperature in other collectors, which may need to be defocused as a security measure, wasting energy and stressing collector components. The problem of defocusing has been widely addressed, for example, using model predictive control (MPC) (Camacho and Bordons, 1995), in works such as (Frejo and Camacho, 2020; Sánchez et al., 2019, 2018b). Authors as Masero et al.; Sánchez et al. (2018a) propose controlling the flow in each collector by using valves, but this requires an estimation of the spatial distribution of the DNI. Since plants typically have few pyrheliometers, other information sources are combined, e.g., images (Radovan and Ban, 2014) and weather information (Zhang et al., 2018). Also, in (Gallego and Camacho, 2012), the spatially distributed irradiance is obtained from the model of the plant using a Kalman filter. Finally, this thesis, which has been performed in the framework of the OCONTSOLAR project,1 considers low-cost direct normal irradiance (DNI) sensors mounted on unmanned aerial and ground vehicles as a mobile robot sensor network (RSN) (Huang et al., 2019; Julian et al., 2012a; Koubâa and Khelil, 2014), a particular case of wireless sensor networks (WSN) (Akyildiz et al., 2002; Aydin et al., 2019) to collect environmental variable in the PTC-CSP (Lilliestam and Pitz-Paal, 2018). Figure 1.1 OCONTSOLAR Project logo. OCONTSOLAR researchers develop new methods to integrate the RSN and its spatially distributed solar irradiance estimation within the control system, improving its predictions to achieve safer and more efficient plant operation (Aguilar-López et al., 2022), as can be seen in Fig. 1.2. In this scheme, a DNI estimator processes the data from pyrheliometers and other sources of information such as cameras, meteorological information, and the measurements of the collector sensors and the robots. This module feeds the PTC-CSP control system, which is detailed in Fig. 1.3, and the robot fleet block shown in Figure 1.4. The former computes the plant states from the combination of the radiation maps with the flow and temperature measurements at the solar collectors and calculates control actions (valves opening, pump power and defocusing of collectors); the latter uses the DNI maps and the DNI certainty map to generate a set of tasks that will be assigned to the robots depending on their state, e.g., position and battery level, and characteristics. The proposed control scheme can be seen in Fig. 1.5. Finally, this thesis focuses on the robotic fleet and its control, and on the use of the measurements gathered to estimate the DNI maps. 1Optimal Control of Thermal Solar Energy Systems H2020 ADG-ERC Project (Grant Agreement 789051) 1.1 Multi-Robot Systems for Monitoring Solar Plants 5 Figure 1.2 OCONTSOLAR scheme. Figure 1.3 PTC-CSP control scheme. Figure 1.4 Robot fleet control scheme. 1.1 Multi-Robot Systems for Monitoring Solar Plants Multi-robot systems (MRS) are groups of robots that collaborate to fulfill a set of tasks (Cortés and Egerstedt, 2017; Rizk et al., 2019) difficult to address using a single robot. In this thesis, we focus on MRS to provide a RSN, for mapping and monitoring environmental variables (Bolstad, 2016). Indeed, MRS have been used lately in many applications. For instance, MRS can perform surveillance tasks (Ahmadi and Stone, 2006; Stump and Michael, 2011), e.g., to patrol marine and urban environments (Mendonça et al., 2016; Semsch et al., 2009); in logistics, cooperative MRS can find and pick up objects in unknown locations (Bähnemann et al., 2017) and maximize the number of packages transported (Farinelli et al., 2017); in precision agriculture (Ribeiro and Conesa-Muñoz, 2021), MRS, inspect the terrain and carry out the corresponding treatments, which can 6Chapter 1. Introduction Figure 1.5 Proposed control scheme. include weed and pest control (Gonzalez-de Santos et al., 2017); MRS can also be used for filming sports events with the aim of maximizing viewer’s satisfaction (Zema et al., 2017), and for inspection, e.g., to monitor the state of the blades of wind turbines (Franko et al., 2020); MRS also build 3 D maps using the data supplied by 2 D laser range finders and combining scan-matching with sample-based probabilistic localization (Nex and Remondino, 2014; Thrun et al., 2000); MRS have also been proposed to manipulate flexible payloads in the context of automotive assembly (Mills and Ing, 1996); Likewise, MRS can provide communications support in areas with weak infrastructures (Wang et al., 2020). Other applications of RSNs that are more relevant in the context of the present thesis are given in spatial estimation (Roldán et al., 2016), where one unmanned ground vehicle (UGV) and one unmanned aerial vehicle (UAV) are used to monitor temperature and humidity in a greenhouse, and in (Conesa-Muñoz et al., 2016), where aerial and ground vehicles are respectively used to gather environmental information and perform interventions. Also, in (Zhang and Leonard, 2010), a cooperative Kalman filter is used for both managing a RSN and estimating the state of a static planar scalar field. 1.1.1 Multi-Robot Task Allocation In this context, a common challenge is the so-called Multi-Robot Task Planning (MRTP), where a set of tasks must be assigned to the robots in the most efficient way. This problem has been broadly studied during the last two decades, as can be seen in (Khamis et al., 2015; P. Gerkey and Mataric, 2003a,b; Seenu et al., 2020; Yan et al., 2013), and addresses the allocation of tasks to robots to optimize one or more performance indexes, such as energy consumption and total time spent among others. 1.1 Multi-Robot Systems for Monitoring Solar Plants 7 MRTP is usually subdivided into two different problems: Multi-Robot Task Decomposition (MRTD), which refers to how the missions are decomposed into single subtasks that can be completed by independent robots, and Multi-Robot Task Allocation (MRTA), which refers to the allocation of the tasks and subtasks to the robots. Some relevant concepts in this regard are: • Decomposability: A task is decomposable if it can be divided into a set of subtasks. Also, a task is multi-decomposable if it admits more than one possible decomposition. Attending to the decomposability, we can classify tasks into i) elemental tasks, or atomic tasks if they cannot be divided; ii) simple tasks, if they can be decomposed into elemental or decomposable simple subtasks; iii) compound tasks, if they can be divided into simple or compound subtasks with the condition that there is only one fixed full decomposition; and iv) complex tasks, if it is composed of subtasks that can be simple, compound or complex. • Constraints: These are limits that restrict the space of feasible solutions, e.g., capability constraints define which robots are capable of performing which tasks, simultaneity and non-overlapping constraints specify whether two or more tasks need to be fulfilled at the same time, precedence constraints specify whether a certain task must be carried out before another one, and path constraints refer to the possible paths a robot can follow. • Utility: It measures the difference between the suitability of assigning task tj to robot ri and the cost of performing the task, which becomes infinite if the robot cannot perform the task. According to the taxonomy proposed in (Gerkey and Matarić, 2004; Gerkey and Mataric, 2004), MRTA problems can be classified deppending on: • The number of tasks each robot can perform at a time: they can be Single-Task robot (ST) or Multi-Task robot (MT) problems. • The number of robots needed to perform a task: they can be Single-Robot (SR) or Multi-Robot task (MR) problems. • The availability of information to plan future allocations: they can be Instantaneous Assignment (IA) or Time-Extended Assignment (TA) problems. This taxonomy has been further developed in (Korsah et al., 2013) by considering the interdependencies between tasks, with no/in-schedule/cross-schedule/complex dependencies (ND/ID/XD/CD) problems, taking into account whether there is no/intra-schedule/interschedule dependency between simple tasks, or inter-schedule dependency between complex tasks. Note that dependency is highly related to the aforementioned utility concept. Also, another classification deals with the differences on the robots used. MRTA problems are homogeneous if all the robots have the same characteristics. Otherwise, they are heterogeneous MRTA problems. Likewise, two types of approaches can be distinguished within MRTA solutions: centralized and decentralized or distributed. While the decentralized approach tends to be 8Chapter 1. Introduction more scalable and robust to communication failures because it divides the overall problem into pieces containing partial information, the centralized one benefits from its access to broader information to generate allocations closer to the optimal one, i.e., the centralized approach improves performance, but it increases the computational load. Furthermore, centralized approaches are less robust against failures and attacks; and are less able to adapt themselves to dynamic environments (Yan et al., 2013). One halfway possibility between centralized and distributed approaches is that of clustering (Chanfreut et al., 2021; Dörfler et al., 2014; Rahwan et al., 2015; Wang et al., 2018), which seeks a trade-off between performance and computation burden. In this way, sets of tasks and robots are partitioned into loosely coupled clusters, which can be efficiently managed in parallel. Performance also depends on the algorithm used. The most commonly used algorithms for these problems are optimization-based and market-based algorithms. On the one hand, optimization-based approaches (Horst et al., 2000) focus on finding the optimal solution of a constrained problem (Nam and Shell, 2014). These algorithms tend to be centralized and can lead to NP-Hard problems as in the Optimal Assignment Problem (OAP) (Gale, 1989; Kuhn, 1955a). In order to avoid high computational costs, some metaheuristic algorithms (Juedes et al., 2004; Yang, 2011) have been used to perform the optimization, as in (Farinelli et al., 2017), where graphical models such as distributed constraint optimization problems (DCOPs) are explored. Within these algorithms there are two main subtypes; namely, deterministic and stochastic: • Among deterministic algorithms, we can find works based on the gradient and the Hessian methods such as (Atay and Bayazit, 2006; Darrah et al., 2005; Xue et al., 2019), where the MRTA is solved using mixed-integer linear programming algorithm (MILP). • By contrast, among stochastic algorithms we can find different types, e.g. Jones et al. (2011a); Jose and Pratihar (2016); Panchu K. et al. (2018); Tolmidis and Petrou (2013) proposed genetic algorithms (GA) for the MRTA problem. Also, authors as Huang et al. (2018); Xue et al. (2020) adapt the well-known cuckoo search to solve it. Finally, particle swarm optimization (PSO) (Asma and Sadok, 2019; Choudhury and Biswal, 2011; Li and Ma, 2008) and ant colony optimization (ACO) (Li et al., 2017; Wang et al., 2012) have been widely studied. On the other hand, market-based algorithms (Dias et al., 2006) are inspired in economy and provide a simple way of allocating resources but imply explicit communication among the agents. These algorithms are more robust and scalable than optimization-based ones. However, the communication cost can be too high. Moreover, since there is no agent that has full information, the final allocation is not guaranteed to be an optimal one. The most widespread ones are auction protocols, where agents indicate their interest in resources by bidding. Indeed, auction algorithms have already been used to solve the MRTP problem (Gerkey and Mataric, 2002), because MRS are cooperative, i.e. robots are not typically programmed to be selfish or lie. Also, these algorithms allow us to take into account the availability of different types of resources and the cost of communications (Lee, 2018) and use consensus to reach an intermediate solution (Choi et al., 2009). Furthermore, auctions algorithms are able to deal with complex tasks (Zlot and Stentz, 2006b). However, 1.2 Outline 9 these algorithms may share some of the problems found in centralized approaches, e.g., sensitivity to failures in communication (e.g., there must be an auctioneer for the bidders). Therefore, in works such as (Gerkey and Mataric, 2002; Turner et al., 2018; Zhao et al., 2016), there is a great emphasis on monitoring the progress of the allocated task and in task reallocation whenever in case of need. This way, any failure in communication can be detected and solved. Other works that use a market-based algorithm avoid the auction mechanism by introducing a market where any robot can exchange tasks with the others (Tang and Parker, 2007; Zlot and Stentz, 2006a). After introducing these concepts, we can properly define the problem that has been addressed in this thesis. In particular, we have considered a homogeneous robot fleet composed of UAVs and UGVs to gather DNI data all over the PTC-CSP, which has a structured layout with known obstacles. This problem can be classified as ST-SR-MRTA because a robot can only take a measurement at a time and only one robot is required to take a measurement. Furthermore, it can also be classified as a XD problem because the utility that a task has for a robot depends not only on its own schedule but on the schedule of the rest of the robots. Finally, in this thesis, a central layer that performs the MRTA has been considered. Thus, the strategies used to approach the problem are centralized. 1.2 Outline This thesis has been divided into four parts. Part I is composed of Chapters 1 and 2, and surveys the literature, formalizing the problem setting. Part II focuses on blocks B6 and B7of Fig. 1.5, i.e., MRTA and robot dynamics. More specifically: • Chapter 3 deals with the MRTA problem (see Fig. 1.5) and presents an IA approach. • Chapter 4 introduces two TA approaches: a branch and bound algorithm (BaB) and a GA for MRTA problems. • Chapter 5 proposes a clustering algorithm based on cooperative game theory to split large MRTA problems into multiple smaller MRTA problems. • Chapter 6 presents a predictive control scheme with a receding horizon strategy for the case of dynamic tasks. Part III focuses on block B1 of Fig. 1.5, i.e., the DNI estimation. The outline of this part can be seen below: • Chapter 7 considers the task generation and the DNI estimation problem, developing an algorithm that accounts for the probability distribution of the DNI at each point of the plant and its update from new measurements using Bayesian inference. • Chapter 8 employs spatio-temporal kriging to perform the DNI estimation from a mesh of sensors. Finally, Part IV contains Chapter 9, where the conclusions are discussed and some future lines of work are listed. 16 Chapter 2. Problem Setting Some of the elements of the plant previously mentioned can be considered as obstacles for robot navigation: manifolds are unsurmontable for UGVs because they are at ground level, and require the setting up of bridges or passing areas; also, the possible fall of an UAV represents a serious risk for the integrity of parabolic mirrors, so they are not allowed to fly over them; finally, the whole area of the buildings is reserverd for human use, so no robots are allowed into it. 2.2.1 Solacor I, El Carpio Many of the chapters are based on a section of the plant Solacor I in El Carpio, Spain (Abengoa), whose layout can be seen in Fig. 2.4 and Fig. 2.5. Figure 2.4 Section of the real Thermosolar power plant of El Carpio (Spain) PTC-CSP (from Google Maps). Figure 2.5 Layout considered. This 30 MW thermosolar plant covers 63 ha ( 1180 ×540 m), containing its operation 2.3 Environment Discretization 17 zone in the upper-center of the plant. Likewise, cold and warm HTF manifolds run parallel from the southern part of the operation zone to the extremes of the plant, dividing the plant into three sectors. There are 18 loops in the south sector of the plant and 8 loops in each of the north sectors, and a bridge for the UGVs to go from the south to the north at each side of the operation zone. Each loop has 4 collectors in “U" shape because loops start and end at the manifolds. 2.3 Environment Discretization The environment of the PTC-CSP that must be covered by the robots is denoted by Q⊂R2 . The environment contains a finite number of obstacles for UAVs, OA=∪OA i=1OA i , and those for UGVs, OG=∪OG i=1OG j , where OA i⊂Q for all i=1,2,...,OA , OG j⊂Q for all j=1,2,...,OG, OA i∩OA i′=/0 ∀i=i′,OG j∩OG j′=/0 ∀j=j′.(2.1) These obstacles are assumed to be 2-D convex polygons. For convenience of the subsequent discussions, we assign IDs to vertices of OA i and OG j . The ID for l -th vertex of OA i ( OG j ) is denoted by rA il ( rG jl ), and its position coordinates are denoted by ρA il ∈Q ( ρG jl ∈Q ). The collection of ID of all vertices for all OA 1,...,OA OA ( OG 1,...,OG OG ) is denoted by VA(VG). Let us finally partition the environment Q into a finite number of cells, namely C , labeled C={c1,c2,...,cc,...,cC}.(2.2) The position coordinates of the center of cell cjare denoted by qj∈Q. Alternatively, these cells can be labeled using a couple of subindexes i=1,2,...,Nr and j=1,2,...,Ncto denote the row and the column of the cell, C={c11,c12,...,ci j,...,cNrNc}.(2.3) 2.4 Tasks In this work we consider a set of T tasks, T={t1,t2,...,tT} , assuming that fulfilling a task, tj , consists of sending any robot to a certain cj and staying there for a given time, ψj . Also, note that the relevance of the task is weighted by ϕj (the larger, the more urgent). Finally, consider that the time to perform a task only depends on the task itself. 18 Chapter 2. Problem Setting 2.5 Non-Euclidean Distance Function In many chapters, we have used distance from robots to tasks as a metric computed using graphs to model the possible robot trajectories (Fentanes et al., 2015). To this end, we define a weighted (undirected) graph Di associated with robot ri over the environment Q . Let us now take a UAV ai ( ri∈RA ). Define the node set NA i as NA i={ri}∪RA∪C , where RA is the set of vertices of the obstacles for UAVs and C is the set of IDs for the cells. Let us also define NG i for a UGV gi ( ri∈RG ) in the same way and the edge set EA i . For this purpose, we define a map X that charts elements of NA i and NG i to their position coordinates, namely X(n) =        pi,if n=ri qj,if n=cj ρA il ,if n=rA il ρG jl ,if n=rG jl .(2.4) Using this map, the edge set EA i is defined as below. Take a pair of any two elements ni and nj(ni=nj) in NA i. Then, {ni,nj}∈EA iiff βX(ni)+(1−β)X(nj)/∈OA∀β∈[0,1],(2.5) which means that the line connecting X(ni) and X(nj) does not intersect with the obstacle for UAVs. The edge set EG ifor a UGV ri∈RGis also defined analogously. The weights wA i j ∈WA ( wG i j ∈WG ) corresponding to an edge {ni,nj}∈EA i ( EG i ) are defined by the Euclidean distance ∥X(ni)−X(nj)∥2 . In summary, an undirected weighted graph Diis defined for each ri∈Vas Di= (NA i,EA i,WA i)(2.6) if ri∈RA, and as Di= (NG i,EG i,WG i)(2.7) if ri∈RG. Let us now define the distance between ri and qj following the definition of (Elgindy and Mitra, 1994; Mitchell, 2017; Wang et al., 2011). Then, the shortest path from ri to cj along with the graph Di can be first computed by using any computationally efficient algorithm and once the path is fixed, the travel distance is analytically formulated as a function of pi (i.e., the position of ri ), and qj and it is denoted by d(pi,qj) . Note that we also denote the distance between two cells by d(qj,qj′). 2.6 Artificial DNI Data Generation 19 2.6 Artificial DNI Data Generation To test the proposed approaches, we need data of DNI throughout the PTC-CSP. However, obtaining such data is very difficult as it requires a very fine mesh of sensors. Thus, these data have been generated artificially. With this purpose, clouds are assumed to be composed of a cluster of ellipsoids in random directions contained in a larger ellipsoid with random dimensions following a Gaussian distribution according to (Kulemin, 2003). An example of these clouds can be seen in Fig. 2.6. Figure 2.6 Example of clouds above the power plant. We consider wind as a vectorial field in the volume contained in the limits of the plant area. Note cumulus clouds usually appear between 500 m and 2000 m. It is important to keep in mind that the higher the altitude the faster the wind becomes, following the well-known Hellmann equation: vW vo W =hW ho Wα (2.8) where vW is the speed of wind at height hW and vo W is the speed of wind at the measurement height ho W . α is the Hellmann coefficient, which depends on the type of terrain, and is assumed to be 0.2here due to the moderately rough terrain. For greater realism, and given the turbulent nature of the wind, we will consider that in each one of the mesh cells there will be a disturbance N(0,σV)∈R3 . Since we have the velocity field, we can interpolate the velocity of each cloud in each time and then use 20 Chapter 2. Problem Setting Spencer equations (Spencer, 1971) to obtain solar rays at the center of each cell and their interference with clouds. Since nominal clear-sky irradiance can be computed with clear sky models, irradiance losses can be modeled by a cloud factor (CF) in the range 0−1 in each measurement spot, creating the CF grid of Fig. 2.7. If CFi j =0 , spot i j receives nominal DNI, i.e., there is no shadow, whereas CFi j =1 implies that spot i j receives 0 DNI, i.e., is totally shadowed. Finally, discretized real values are denoted as CFR i j. 1 … j ... Nc Nr i 1 CFR Nr1 CFR Nrj CFR NrNc CFR i1 CFR ij CFR iNc … … … … … … CFR 11 CFR 1j CFR 1Nc Figure 2.7 CFRgrid. By generating random clouds and various wind fields as described above, we generate different data sequences as shown here.3 2.7 Objectives To carry out the optimization, the proposed cost function must take into account multiple criteria: • Time employed to complete each task, i.e., the amount of time elapsed between the allocation and the completion of a task. •Total distance traveled by each robot. •Feasibility of the allocation from an energetic viewpoint. • Each task may have a different weight as not all the tasks are equally urgent. Moreover, it may be more desirable to use a certain robot than another, which means that the distance traveled by the robots also has different weights. 3Video in https://youtu.be/LkV-yv38ONc 2.8 Cloud Estimation 21 • It is also necessary to take into account that allocations must be feasible, i.e., a robot must not only have enough energy to accomplish its allocated tasks but also to go to a charge station after completing them. 2.8 Cloud Estimation Sensors can be fixed or mounted onto mobile robots that take DNI measurements and have a predefined resolution and range, i.e., they can only take values from an alphabet A . We denote the set of measurements as M={m1,m2,···,mM} . Note that each measurement is a set containing the cell where they were taken, cj(mi)∈C , the time when they were taken, t(mi)∈[to,tf] , being to the start time of the simulation and tf the end time of the simulation, and the measured value, u(mi)∈A , i.e., mi={cj(mi),t(mi),u(mi)} . The goal is to obtain CFEst using set Mso that E= Nr ∑ i=1 Nc ∑ j=1 tf ∑ t=to|CFEst i jt −CFR i jt |(2.9) is minimized. To do so, we will equivalently minimize Et=1 Nr·Nc Nr ∑ i=1 Nc ∑ j=1|CFEst i jt −CFR i jt |,t∈[to,tf].(2.10) Then, our goal is to obtain CFEst i j(t=t⋆) (with t⋆ the current time instant) from the previous measurements, i.e., those that fulfill t⋆−t(mi)≥0 . Furthermore, forecasting is equivalent to obtaining a certain CFEst i j(t=t⋆) using only measurements that fulfill t⋆−t(mi)≥F , with Fthe forecasting time. Finally, we can use the label of c∈Cas an alternative to write (2.9) and (2.10) as E= tf ∑ t=to Et(2.11) Et=1 Nr·Nc C ∑ c=1|CFEst ct −CFR ct |,t∈[to,tf].(2.12) 2.9 Further Comments on the Problem Setting The robot fleet described in 2.1 will be the same for all chapters except Chapter 8. The parameters listed in Table 2.1 are not used in Chapter 5, since in this chapter different random parameters are considered for all the robots. Section 2.5 is relevant for all the chapters that consider the PTC-CSP plant layout described in Section 2.2, i.e., Chapter 3, Chapter 4 and Chapter 7. Chapter 5 and Chapter 6 22 Chapter 2. Problem Setting consider environments with no obstacles for simplicity. However, disregarding the presence of obstacles, the plant has been discretized in all the chapters as explained in Section 2.3. Finally, tasks and objectives introduced in Section 2.4 and Section 2.7 are used in Chapter 3, to Chapter 7. Likewise, artificial cloud data generation and cloud estimation objectives described in Section 2.6 and Section 2.8 respectively are used in Chapter 7 and Chapter 8. Part II Multi-Robots Task Allocation 23 3 Event-MILP Algorithm for IA MRTA Problem in the PTC-CSP This chapter is taken from (Martín et al., 2021), where an event-MILP-based algorithm is proposed to solve the task allocation problem of a RSN. A simulation environment taking into account the behavior of the robots and the structure of the solar plant was built to test the proposed algorithm, which was then compared with traditional methods such as OAP. The rest of the chapter is organized as follows. In Section 3.1 the problem under consideration and the objectives pursued are presented. Besides, the power solar plant model and the allowed paths for each type of vehicle are introduced. The task management algorithm is introduced in Section 3.2. The simulation environment including the vehicle models and their behavior are shown in Section 3.3. Results from simulations are analyzed in Section 3.4 and finally, some conclusions are drawn in Section 3.5. 3.1 Case-Study As introduced in Chapter 2, we have a set of robots R=RG∪RA that must operate continuously in a PTC-CSP that has been discretized in C , accomplishing tasks to be introduced in batches, B , periodically. Furthermore, robots must operate avoiding the corresponding obstacles in OG and OA . With this purpose, routes and distances have been computed as was described in Chapter 2.5. Here, the set of tasks T changes constantly. Tasks that have been already allocated to a robot are removed from the set, while every time there is a new batch of tasks we include them in the tasks set, T′=T∪B . Then, the T tasks must be allocated to the robots periodically and the goal of the fleet of robots is to accomplish the tasks in the most efficient way possible attending to the criteria established in Chapter 2.7. The PTC-CSP considered is the one described in the layout of Chapter 2.2. Note that, from the UGV perspective, the arrangement of obstacles leads to the division of the field into three sectors (marked with orange boxes) separated by two bottlenecks, which can be observed in Fig. 3.1. Also, the C=434 central positions ( qj ) in which the solar plant has 25 32 Chapter 3. Event-MILP Algorithm for IA MRTA Problem in the PTC-CSP In case all charge stations are being used, this robot will remain “Stopped” until a charge station is free again. 10. In the case there are no more tasks for UAVs they land and in the case they are currently on land when a task is loaded, the first thing that they must do is to reach their reference altitude. The simulation environment can be observed in a more schematized way in the flowchart in Figure 3.4 Figure 3.4 Simulation environment scheme. The robot models used are simplified models to obtain a lighter and faster simulation. It is supposed that each robot has a perfect low level controller so that they can be considered 3.4 Results 33 as points that move on the space. The landing and take-off are also considered for UAVs, landing when they are not performing any task and taking off when a new task is allocated to them. Also, energy consumption and recharging are taken into account. The parameters considered for UGVs and UAVs in the simulations are depicted in Chapter 2.1. The simulation environment depicted here has been programmed in Matlab ® and intlinprog function from the Optimization Toolbox has been used in order to solve the cost function (3.1). 3.4 Results The proposed algorithm has been tested considering 2 UGVs and 3 UAVs and a duration of (tf−to) = 12 hours that corresponds to an average working day on a thermosolar plant. During these 12 hours a total of 420 measurement tasks appear. These measurement tasks have been generated randomly with the same seed so that they are the same for the compared algorithms. The computation time required to perform these simulations has been 4 hours and 5 minutes for the OAP based approach and 4 hours and 15 minutes for the proposed approach. The weights for the proposed algorithm have been tuned fixing the values for the weights related to the cost of allocating a robot i to a task j , F1 and F2 , and varying the weights related to not allocating task j , i.e., F3 and F4 , equations (3.6) and (3.7) respectively. The proposed algorithm has been compared with the case in which not allocating a task is not allowed unless there are not enough robots to allocate all the tasks, i.e., F3=∞ as stated in Section 3.2. All the values for the weights are in Table 3.1. Table 3.1 Weights. Values can differ for UGVs and UAVs. Proposed algorithm OAP algorithm F1UGV: 1/ UAV: 0.05 F2UGV: 1.1−0.1·bi/ UAV: 1.1−0.1·bi F30.05 ∞ F40.15 - To analyze the performance, several indexes usable for any MRTA algorithm have been defined. They consider different points of view such as distance, energy, dispersion, etc, and their values are represented in Tables 3.2 to 3.9 (the best results of each index is highlighted in bold. In case neither of them is highlighted in bold it means that any of them can be considered better than the other one). •The distance traveled by the robots, since distance is one of the parameters used –In1is the total distance traveled by all the robots. –In2 is the sum of the remaining distances that each robot still needs to travel to complete all the tasks stacked in its queue. 34 Chapter 3. Event-MILP Algorithm for IA MRTA Problem in the PTC-CSP –In3 is related to the total distance that the robots travel, i.e., the sum of the remaining distances that each robot still needs to travel to complete all the tasks stacked in its queue and the distance that they have already traveled. –In4is the percentage of distance traveled by UGVs. –In5is the percentage of distance traveled by UAVs. Table 3.2 Indexes relate to distance . Proposed algorithm OAP algorithm In1250204.02 245949.86 In21305 5109 In3251509.02 251058.86 In417.45 18.20 In582.55 81.80 •The dispersion of the robots around the map which can be calculated by: disp =1 R·(tf−to) R ∑ i=1 tf ∑ t=to∥¯pi−pi(t)∥,(3.9) where ¯pi is the average position of ri among the complete simulation. These are the most interesting indexes since the main aim of leaving tasks unallocated is having the possibility of allocating them in later batches, performing better allocations 1. In6is related to the dispersion of all the robots. 2. In7is related to the dispersion of the UGVs. 3. In8is related to the dispersion of the UAVs. Table 3.3 Indexes relate to dispersion. Proposed algorithm OAP algorithm In6406.44 442.05 In7311.58 296.95 In8469.68 538.78 •The energy, which is related to the traveled distance. 1. In9 is the mean of the energy of the robots. Note that the less this index value is, the more energy the robots will have used. 2. In10 is related to the total energy that robots have at the end of the simulation. 3.4 Results 35 Table 3.4 Indexes relate to energy. Proposed algorithm OAP algorithm In959.13 58.25 In10 33.50 49.58 •The stacked tasks during the simulation. 1. In11 is the maximum number of stacked tasks during the simulation. 2. In12 is the mean number of stacked tasks during the simulation. Table 3.5 Indexes relate to stacked tasks. Proposed algorithm OAP algorithm In11 25 25 In12 3.83 4.21 •The number of tasks that have been performed during the simulation. 1. In13 is the number of tasks of any type completed by all the robots. 2. In14 is the number of charge tasks completed by all the robots. 3. In15 is the number of charge tasks completed by UGVs. 4. In16 is the number of charge tasks completed by UAVs. 5. In17 is the number of measurement tasks completed by all the robots. 6. In18 is the number of measurement tasks completed by UGVs. 7. In19 is the number of measurement tasks completed by UAVs. Table 3.6 Indexes relate to the number of tasks performed. Proposed algorithm OAP algorithm In13 441 432 In14 25 24 In15 0 0 In16 25 24 In17 416 408 In18 78 84 In19 338 324 •The number of non allocated tasks. 1. In20 is the mean number of unallocated tasks during the simulation. 2. In21 is the number of unallocated tasks at the end of the simulation. 36 Chapter 3. Event-MILP Algorithm for IA MRTA Problem in the PTC-CSP Table 3.7 Indexes relate to the non allocated tasks. Proposed algorithm OAP algorithm In20 1.1179583 0.0009726 In21 2 0 •The time that was employed to serve the tasks. 1. In22 is the mean time employed to serve the tasks, taking into account both types of tasks. 2. In23 is the mean time employed to serve the charge tasks. 3. In24 is the mean time employed to serve the measurement tasks. Table 3.8 Indexes relate to the time to serve tasks . Proposed algorithm OAP algorithm In22 250.51 252.92 In23 2396.59 2419.56 In24 121.54 125.47 • The time employed to solve the allocation problem, which is one of the main parameters while evaluating the general performance of the algorithm. 1. In25 is the mean time employed to solve the allocation algorithm. 2. In26 is the maximum time employed to solve the allocation algorithm. 3. In27 is the total time employed to solve the allocation algorithm during the simulation. Table 3.9 Indexes relate to the solving allocation times. Proposed algorithm OAP algorithm In25 0.01 0.01 In26 0.22 0.21 In27 6.29 3.92 It is important to emphasize that although the tasks are the same for all cases, the allocation may affect the order in which they are carried out. Note as well that, even though the best values for each considered index have been marked in bold in Tables 3.2 to 3.9, the indexes are coupled. Thus, it is important to analyze the values of the indexes using a global perspective. For example, it could be misleading to compare the indexes related to the distance traveled without taking into account the number of tasks performed and other variables. In Fig. 3.5 an analysis of the results in Table 3.10 is represented considering that the indexes are the relativized values with respect to the OAP. 3.4 Results 37 Table 3.10 Comparison of the results. Proposed algorithm OAP algorithm Proposed / OAP 1: Total distance / total tasks (In3/In13)570.32 581.15 0.9813 2: Stacked tasks / total tasks (In11/In13)0.05669 0.05787 0.9796 3: Mean energy consumed / total tasks ((100 −In9)/In13)0.09268 0.09664 0.9944 4: Mean energy consumed / total distance (((100 −In9)/In3)) 0.0001625 0.0001663 0.9771 5: Measurement tasks / charge tasks (In17/In14)16.64 17 0.9788 6: Total dispersion / total tasks (In6/In13)0.9216 1.023 0.9007 7: Mean time to serve measurement tasks (In24)121.54 125.47 0.9687 8: Mean time to serve charge tasks (In23)2396.59 2419.56 0.9905 9: Total time allocating tasks / total tasks (In27/In13)0.0142 0.0091 1.5718 (a) General view. (b) In detail. Figure 3.5 Indexes Radar Chart of the algorithm proposed and OAP. Based on the results of the different indexes, and taking as a reference the OAP algorithm, the following conclusions can be obtained: • Generally, using the proposed approach made the dispersion of the robots decrease. This was in part as a result of a better grouping of the tasks, which was the main aim of the proposed algorithm. However, in some cases, due to a bad tuning of the weights that affect ϑ , the time that the algorithm let some tasks unallocated can have affected badly the mean time in which measurement tasks were carried out, or even the number of tasks performed. • In general, the lower the ϑ are, the more the indexes related to time employed to solve the allocation increase. This is logical since these times depend mostly on the number of variables, their associated costs and the number of auctions carried out. 38 Chapter 3. Event-MILP Algorithm for IA MRTA Problem in the PTC-CSP • It has to be taken into account that the fact that tasks are reordered in the queue of each robot every time a new task is allocated to them improves the performance of those cases in which the algorithm favors allocating tasks ( ϑ are high from the beginning and increase fast) including the OAP. • When the costs of not allocating tasks are lower, the mean number of unallocated tasks, In20 , increases. Besides, the maximum number of stacked tasks, In11 , is lower as well. • By tuning the weights F3 and F4 properly, not only has the dispersion of the robots been minimized but the mean time on serving charge and measurement tasks has also been reduced 0.95% and 3.13% respectively. It can also be seen that with this tuning, more measurement tasks than in the OAP have been fulfilled and that the ratio of distance traveled per measurement task is bigger as well. • Every time that the pile of stacked tasks increases and that UGVs fulfill tasks, it is as a consequence of having the UAVs in “charge mode”. This can be seen in Figure 3.6 and in Figure 3.7, since UGVs are significantly slower than UAVs and since the algorithm favors allocating tasks to UAVs. • Depending on what goals are being pursued, the weights can be tuned to prioritize different objectives such as using less energy, using a certain type of robot better or fulfilling the tasks faster. In summary, although the algorithm proposed requires more computation time to carry out an allocation since it deals with a larger amount of alternatives (tasks can remain unallocated and unallocated tasks add up to general tasks for future allocations), it can be concluded that it achieves its goal of grouping the tasks by zones, thus achieving a better performance in the course of the full simulation. Figure 3.6 SOC and stacked tasks of the robots in the proposed algorithm. 3.5 Conclusions 39 Figure 3.7 SOC and stacked tasks of the robots in the OAP algorithm. 3.5 Conclusions The algorithm proposed in this chapter allocates tasks in an asynchronous way, considering the option of not allocating tasks with an associated cost. The cost of not allocating a task varies with the time that those tasks have been unallocated. The utilities that have been used by each robot to measure its cost to perform a certain task depend on the distance that it has to travel, the distance that it will be from a charge station once the task is completed, and the energy level that it has at the moment of receiving the allocation. This algorithm has been tested in simulations, considering a simulation time of twelve hours and using several indexes that have been proposed to compare any MRTA allocation problem. It has been found out that, as long as the weights are well-tuned, having the possibility of not allocating tasks can involve benefits in the grouping of tasks so that they are performed in a more efficient way. However, a badly chosen cost on not allocating tasks can be detrimental as tasks may not be allocated even in a case where performing them could be better than not doing so. Although the proposed algorithm means an increase in the allocation solving time with respect to the OAP, it may improve the efficiency in terms of energy consumption and distance traveled and the mean time for serving the tasks. 4 BaB and GA for TA MRTA Problem in the PTC-CSP In Chapter 3 we tried to improve the overall performance of the robot fleet over time, but it is difficult due to the limitations of IA MRTA. The results of this chapter, published in (Martin et al., 2021a), deal with the TA MRTA problem. Thus, we take into account the route that the robots will follow when solving the allocation considering the criteria established in Chapter 2.7. Also, the computational burden of the problem is drastically reduced by using discrete variables, so that obtaining an allocation can be addressed by a BaB in small problems and by a customized GA for larger problems. The design of the GA crossover and mutation strategies ensures that the descendant allocations of each generation will maintain a certain level of feasibility, reducing greatly the range of possible descendants, and accelerating their convergence to a suboptimal allocation. The proposed MRTA algorithms are simulated and analyzed in the context of the PTCCSP of Chapter 3, for which the spatially distributed DNI is estimated using a heterogeneous fleet of both UGVs and UAVs. The performance and response of the proposed algorithms are assessed in for four different scenarios. The results show that the GA provides nearly optimal solutions in much less computation time for large problem instances. Likewise, the trade-off between computation time and accuracy can be set by tuning parameters of the GA. This chapter is organized as follows. The variables used to formulated the problem are described in Section 4.1 and in Section 4.2 the cost function is presented. In Section 4.3 and Section 4.4, the branch and bound and the genetic algorithm developed to approach the problem are presented. In Section 4.6 the problem is particularized for a case of study based on a thermosolar power plant and a fleet of robots composed of UAVs and UGVs. The results of proposing several scenarios for the optimization of the problem are analyzed in Section 4.7. Finally, the conclusions and future works are drawn in Section 4.8. 41 48 Chapter 4. BaB and GA for TA MRTA Problem in the PTC-CSP a set boundary, until the optimal solution is found. BaB has been used in some works to address the MRTA problem, as in (Kartal et al., 2016), where it is used in combination with a Montecarlo search tree (Browne et al., 2012). However, it has always been applied to the MRTA problem using binary variables for the task assignments while in this work it is applied using discrete variables (which allows to substantially decrease the number of leaves on the search tree). More specifically, in the present work, the use of the BaB is proposed for solving smallscale (less than 105 solutions according to (4.2) ) MRTA problems formulated as in (4.3) . The algorithm generates all the possible feasible combinations for each k -th allocated task removing the ones that show a cost function higher than the current upper bound for the optimal cost function. For the remaining task allocations the corresponding possible combinations for the following (k+1)th allocated tasks are generated comparing their new cost function with an upper bound. The algorithm continues until all the possible optimal task allocations (for any order) have been explored. More specifically, the algorithm applied in this paper is composed of the following steps: 1. Choose an initial estimation of the BaB threshold (J⋆). 2. Then, all the possible task allocations for the firstly allocated tasks for each robot are generated, i.e., the new “branch" is done by adding a single task to the task queue of each robot (in every possible way). A simple example with R=2 robots and T=2 tasks can be seen in Figure 4.3. The “branch" is done by adding a single task to the task queue of each robot (in every possible way). The entire set of possible allocations for the firstly allocated tasks (composed of 6leaves in this case) is generated. γ11 = 1 γ21 = 0 γ11 = 0 γ21 = 1 γ11 = 2 γ21 = 0 γ11 = 0 γ21 = 2 γ11 = 2 γ21 = 1 γ11 = 1 γ21 = 2 Figure 4.3 Initial allocation tree with R=2robots and T=2tasks. 3. The partial cost function of the generated allocations is computed. This function corresponds to the distances and times of the allocated tasks taking into account that many tasks can be unallocated at this step. The partial cost function used is similar to (4.3) but redefining Ω2(γ) in order to eliminate the penalization for unallocated tasks: 4.3 Branch and Bound Algorithm 49 min γJ⋆(γ) = T ∑ i=1 λi·ρi(γ)+ T ∑ j=1 ϕj·τj(γ)+Ω(γ)(4.13) Ω(γ) = z1·Ω1(γ) s.t.γik ∈T∪{0}∀ i,k In contrast with the original cost function defined in (4.3) , Ω2(γ) is not used since, in this case, also partial assignments which do not fulfill all the tasks are evaluated. Therefore, this partial cost function J⋆(γ) will be lower than an original cost function J(γ) if some tasks are not allocated or are equal to the original cost function if all the tasks are allocated (J⋆(γ)≤J(γ)). 4. The task allocations with a cost function higher than the BaB threshold are discarded as it can be seen in Fig. 4.4, where, assuming that the partial cost function of the first and fourth leaves are higher than J⋆(γ) in the tree shown in Fig. 4.3, these leaves are discarded. γ11 = 1 γ21 = 0 γ11 = 0 γ21 = 1 γ11 = 2 γ21 = 0 γ11 = 0 γ21 = 2 γ11 = 2 γ21 = 1 γ11 = 1 γ21 = 2 Figure 4.4 Bound with J⋆. 5. If any allocation is complete (i.e. all the tasks are allocated) and the corresponding cost function is lower than the current value of the BaB threshold, this value is updated (also discarding the allocations lower than the new cost function) as it can be seen in Fig. 4.5. 6. The allocations which have not been discarded are used to generate a new set of possible allocations including the k -th allocated task for those solutions which have still tasks to allocate, as it can be seen in Fig. 4.6. In this example it is assumed that the cost function of the allocation on the fifth leaf is lower than J⋆(γ)and, as a consequence, the value of J⋆(γ) is updated. Subsequently, the sixth leaf (with a cost function higher than the updated value of J⋆(γ) ) is removed accordingly. Then, the tree is extended considering the task assigned in the second place for the search tree as it can be seen in Fig. 4.6. 7. The algorithm ends if the remaining leaves cannot be increased any more (i.e. if all the tasks are allocated). Otherwise, it returns to step 3 . In Fig. 4.7 the cost function 50 Chapter 4. BaB and GA for TA MRTA Problem in the PTC-CSP γ11 = 1 γ21 = 0 γ11 = 0 γ21 = 1 γ11 = 2 γ21 = 0 γ11 = 0 γ21 = 2 γ11 = 2 γ21 = 1 γ11 = 1 γ21 = 2 Figure 4.5 Cost function is updated. γ11 = 0 γ21 = 1 γ11 = 2 γ21 = 0 γ11 = 2 γ21 = 1 γ12 = 0 γ22 = 2 γ12 = 1 γ22 = 0 Figure 4.6 Extended allocation tree. is updated once more for the tree shown in Fig. 4.6, assuming that the cost function of the complete allocation on the first leaf is lower than the current value of J⋆(γ) (the one corresponding to the last leaf) and, as a consequence, the second and third leaves (with a cost function higher than the updated value of J⋆(γ) ) are also removed. Since only one candidate leaf remains ( γ= [0,1,0,2] ) and this leaf cannot be further increased, the remaining allocation is the optimal one. γ11 = 0 γ21 = 1 γ11 = 2 γ21 = 0 γ11 = 2 γ21 = 1 γ12 = 0 γ22 = 2 γ12 = 1 γ22 = 0 Figure 4.7 New cost function update. 4.4 Genetic Algorithm 51 4.4 Genetic Algorithm The BaB presented above exploits computationally for medium and large size problems. As a result, another faster but suboptimal algorithm has to be used for these cases. The GA is a metaheuristic algorithm inspired by the process of natural selection and evolution which relies on crossover, mutations and selection (Goldberg, 2006; Mitchell, 1998). The GA has been applied to similar MRTA problems in works such as (Jones et al., 2011b), where it is used in combination with intra-path constraints in order to solve the allocation of a heterogeneous fleet of robots in a disaster scenario; (Tolmidis and Petrou, 2013), where a multi-objective cost function including both energy and time criteria is used; (Panchu K. et al., 2018), where a multi-objective cost function is also used in a multi-robot task context including time constraints and discrete variables are used to model the problem; and (Liu and Kroll, 2016) and (Jose and Pratihar, 2016), where similar structured environments with discrete positions of the robots and discrete variables are used. In this work, we have designed a GA including the use of specially developed crossover and mutation strategies that generate ordered solutions for the next generations. These strategies generate new allocations that do not introduce unfeasibility from the viewpoint of tasks (allocations that do not fulfill all the tasks or that fulfill a task more than once cannot be fathered by allocations that do) and thus, reduce drastically the number of possible allocations and the computation load. With the same criterion, it has been taken into account that any allocation for a robot, γi= [γi1,γi2,···,γik]⊂γ , including intermediate zeroes is equivalent to the same allocation removing these zeroes and adding them at the end. The crossover and mutation methods presented in this section are novel to the best of our knowledge. An example of the individuals or chromosomes which are used in this paper can be seen in Table 4.3, where the k -th allocated tasks for each robot ri are concatenated. Note that each individual or chromosome corresponds to a full candidate solution to the MRTA problem (i.e., all the tasks assigned to all of the robots). Table 4.3 Gen example. Robot r1··· ri··· rRrobot Gen example 13··· 5··· 0··· 20··· 0··· 0··· 40··· 0··· 0 k12··· k··· T··· 12··· k··· T··· 12··· k··· T Note that, even though the crossover and mutation strategies developed do not generate unfeasible allocations from the viewpoint of tasks, they can still produce unfeasible allocations from an energetic point of view. The GA applied for the considered MRTA problem is composed of the following steps: 1. An initial population for the GA can be pre-generated. 2. The cost function is computed for the the entire set of individuals of the current generation. The fittest individuals are selected as the elite population, which survive automatically into the next generation. 52 Chapter 4. BaB and GA for TA MRTA Problem in the PTC-CSP 3. The algorithm generates a new set of individuals, which are distributed as follows: a) The elite population fraction, FE , defines the percentage of the new population that is considered to be elite individuals. b) The Crossover fraction, FC , defines the percentage of the non-elite new population which is composed of individuals generated by crossover. An example can be seen in Fig. 4.8. Each crossover strategy is composed of: 1 3 2 0 1 0 0 0 0 0 3 4 4 0 2 0 Individual 1 Individual 2 1 3 2 0 1 0 0 0 0 0 32 0 0 0 0 1 4 3 2 4 1 0 0 0 0 3 2 0 0 0 0 New Individual 1 New Individual 2 Robot 1 Robot 2 Robot 1 Robot 2 3.(b).iii 3.(b).iii 3.(b).iv 3.(b).v (Robot 1, 1st order) 3.(b).v (Robot 2, 2nd order) 3.(b).vi Figure 4.8 Example of the crossover strategy with R=2robots and T=4tasks. i. Two elite individuals are randomly selected. ii. One task is randomly selected. iii. The selected task is removed for both individuals. iv. The task allocations are shifted in order to correct the individual obtained in the previous step. v. The selected task is reassigned to each individual but using the allocation considered by the other individual (i.e. for the first individual, the selected task is now allocated for the position originally given by the second individual). vi. The task allocations are shifted again in order to correct the individual obtained in the previous step. c) FM1 is the fraction of individuals generated by mutations, which are generated by using the first kind of mutation proposed. An example can be seen in Fig. 4.9. Each mutation for the first strategy is composed of the following steps: i. One elite individual is randomly selected. ii. One task is randomly selected. iii. The selected task is removed. iv. The task allocation is shifted in order to correct the individual obtained in the previous step. 4.4 Genetic Algorithm 53 v. A new robot and position is randomly selected for the selected task. vi. The selected task is included in the randomly selected position. vii. The task allocation is shifted again in order to correct the individual obtained in the previous step. 1 3 2 00 04 0 Individual 3 1 2 00 04 0 Robot 1 Robot 2 3 2 0 00 04 0 3.(c).iii 3.(c).iv 3.(b).v (Robot 2, 2nd order) 3.(c).iv 3.(c).vi Figure 4.9 Example of the first kind of mutation strategy with R=2 robots and T=4 tasks. d) FM2 is the fraction of individuals generated by mutations, which are generated using the second kind of mutation proposed. An example can be seen in Fig. 4.10. Note that FM2 =1−FM1 by definition, since there are only 2 types of mutation. Each mutation for the second strategy is composed of the following steps: i. One elite individual is randomly selected. ii. Two tasks are randomly selected. iii. The selected tasks are exchanged. 1 3 2 00 04 0 Individual 4 3 2 00 01 0 New Individual Robot 1 Robot 2 3.(d).iii 3.(d).iii Figure 4.10 Example of the second kind of mutation strategy with R=2 robots and T=4 tasks. e) The algorithms end if the number of generations is equal to the maximum number of generations or if the solution has not been improved during a maximum number of stall generations. Otherwise, the algorithm goes back to Step 2). 54 Chapter 4. BaB and GA for TA MRTA Problem in the PTC-CSP For instance, let us take a case in which we have a population size, PS , of 100 individuals. If FE=0.1 To create the next generation we will copy the 10 best individuals from the current generation. Then, considering FC=0.8 , 80% of the 90 remaining individuals, i.e. 72 individuals, must be created by crossover. To do this we take random selections of elite individuals and do crossovers between them until we reach the 72 individuals we need. The remaining individuals 100 −10 −72 =18 are created by mutations. Assuming FM1 =0.5 and FM2 =0.5 half of them, 9 , will be generated by selecting a random elite individual and applying mutation type 1 and the other half by doing the same but applying mutation type 2. Note that the middle zeroes in the gen do not affect the performance and speed of this approach as the algorithm will not produce not potentially optimal solutions. 4.5 Initial Solutions For both algorithms above, it would be beneficial to have a set of initial solutions I . In the case of BaB the best (lowest) cost function obtained within I can be selected as the initial threshold ( J⋆ ), and in the case of GA, I can be used as the initial population. If the desired initial population is smaller than |I| , we can choose the best individuals among I and, if it is larger, the rest of the population is filled with randomly generated potentially optimal solutions. To do this, we permute the set T randomly and then go through all the tasks, selecting a random robot from R and allocating the task to the robot. This way, the resulting individual will not activate Ω2in (4.3). This section presents the set of suboptimal initial solutions that has been used in this work for the application of the BaB and GAs. More specifically, the set of initial solutions is composed of the following subsets: 1. The R initial solutions resulting from solving the shortest path problem for every robot ri: min γi ρi(γi)(4.14) s.t.Ω2(γi) = 0 γik ∈T∪{0},k∈K Notice that this problem is the well-known traveling salesman problem (TSP). TSP is a NP-Hard problem that can be solved in an exact way using Held-Karp algorithm (Bellman, 1961) in a time O(T22T) . However, the nearest neighbor (NN) algorithm (Johnson and McGeoch, 1997), a heuristic greedy algorithm consisting of choosing the nearest non-visited node in each step, is used to generate feasible solutions in a very reasonable time. 2. The R initial solutions result from solving the shortest time problem for every robot ri . The distances considered in these problems are the ones resulting from adding the 4.5 Initial Solutions 55 estimated traveling times between robots and tasks and the operation times required in the corresponding tasks: min γi τi(γi)(4.15) s.t.Ω2(γi) = 0 uik ∈T∪{0},k∈K This problem is addressed in a similar way to the shortest path problem. However, since the criteria are different, the solution of this problem is both feasible and different from the problem mentioned in (4.14). 3. The initial solution resulting from solving the assignment problem using only distances: min γ′ R ∑ i=1 λi·ρi(γ′)(4.16) s.t.Ω3(γ′) = 0 γi(1)∈T∪{0},ri∈R where γ′= [γ11,γ21,···,γR1]⊂γ ,i.e., the discrete variables corresponding only to the first instant, and Ω3(γ′) = 0 is the set of constraints ensuring that all the tasks are carried out (if T≤R ) or that all the robots are assigned to one task (if T>R ). In case there are more tasks than robots ( T>R ) the problem is iterated by considering the allocation in the subsequent kth orders of the robots. Each iteration takes into account the previous allocation of the robots, the remaining tasks and the distance among them and it is solved using Kuhn-Munkres algorithm (Kuhn, 1955b) in a time O(n3). 4. The initial solution resulting from solving the assignment problem using only times: min γ′ T ∑ j=1 ϕj·τj(γ′)(4.17) s.t.Ω3(γ′) = 0 γi1∈T∪{0},ri∈R Note that this problem can be hard to solve directly since the time in which tasks are performed will depend on the assignment in the previous iteration. Thus, in this set of initial solutions, we address the problem similarly to the distance assignment problem but considering the time of reaching the task and performing it instead of the distance. 56 Chapter 4. BaB and GA for TA MRTA Problem in the PTC-CSP All the initial solutions considered can be fast and easily computed by using well-known algorithms. It has to be pointed out that the used set of initial solutions is independent of the subsequently applied control algorithms (BaB and GA). Therefore, the set used can be adapted for each problem. The only necessary requirement is that the initial solutions have to be computed fast enough to allow subsequent computations within the BaB or GA. In order to have a good performance, it is also important that the initial solutions corresponds to the different criteria of the cost function. For example, if a cost function minimizes time and distance simultaneously the algorithms should be fed with initial solutions reflecting an opposing behavior (i.e. minimizing only time or space). 4.6 Case Study The parameters of the robots in this work were given in Table 2.1. For the specific layout, we have used the same PTC-CSP considered in Chapter 3. 4.6.1 Scenarios To test both, the GA and the BaB approaches, a set of different sized scenarios are simulated: • Scenario I (S-I): In this scenario there are R=4 robots ( RG=2 UGVs and RA=2 UAVs) and a set of T=6 tasks, which according to (4.2) correspond to 60480 possible solutions. • Scenario II (S-II): In this scenario there are R=6 robots ( RG=3 UGVs and RA=3 UAVs) and a set of T=8 tasks, which according to (4.2) correspond to 51891840 possible solutions. • Scenario III (S-III): In this scenario there are R=10 robots ( RG=6 UGVs and RA=4 UAVs) and T=9 tasks, which according to (4.2) correspond to 1.7643·1010 possible solutions. In this scenario it has been considered that all the robots are concentrated in an area of the plant and all tasks are concentrated in the opposite area of the plant. • Scenario IV (S-IV): In this scenario there are R=5 robots, ( RG=2 UGVs and RA= 3 UAVs), and T=15 tasks, which according to (4.2) correspond to 5.0685 ·1015 possible solutions. This scenario models a common situation where there are several tasks per robot. These scenarios are determined by the amount of robots and their types, locations, starting energy, delay on being prepared in the starting point, and penalty for using each one; and the amount of tasks, their locations and the operating time that it takes to perform them. The parameters that define these scenarios can be seen in Table 4.4 for robots and in Table 4.5 for Tasks. Note that in Table 4.4s. 4.6.2 Genetic Algorithm tuning The GA has been tuned with the parameters shown in Table 4.6. These parameters have been chosen heuristically. Given the heuristic nature of the GA, it has been run 10 times in each of the scenarios described in the previous subsection. 4.7 Results 57 Table 4.4 Robot parameters in the four scenarios. Type of robot Starting spot (coordinate x, coordinate y)λibi(%) τo(s) S-I S-II S-III S-IV S-I S-II S-III S-IV S-I S-II S-III S-IV S-I S-II S-III S-IV S-I S-II S-III S-IV r1∈RG∈RG∈RG∈RA(227.5, 386) (161.5, 429) (326.5, 176) (1052.5, 48) 2 2 2 1 98 84 32 18 15 4 10 9 r2∈RG∈RG∈RG∈RG(260.5, 176) (821.5, 48) (182.5, 472) (986.5, 386) 2 2 1 1 94 44 92 8 10 0 2 10.5 r3∈RA∈RG∈RA∈RA(293.5, 48) (1100, 412) (62.5, 344) (953.5, 218) 1 1 2 1 30 76 26 53 1 2 3 11 r4∈RA∈RA∈RA∈RA(1151.5, 48) (755.5, 344) (29.5, 302) (359.5, 5) 1 1 1 1 50 40 28 61 0 4 15 5 r5-∈RA∈RG∈RG- (260.5, 525) (392.5, 472) (1151.5, 134) - 1 2 1 - 25 45 76 - 7 2 8.5 r6-∈RA∈RG- - (953.5, 472) (227.5, 302) - - 1 1 - - 88 32 - - 0 2.5 - r7- - ∈RA- - - (491.5, 525) - - - 1 - - - 25 - - - 5 - r8- - ∈RA- - - (161.5, 386) - - - 2 - - - 26 - - - 0 - r9- - ∈RG- - - (29.5, 429) - - - 2 - - - 12 - - - 0 - r10 - - ∈RG- - - (194.5, 472) - - - 1 - - - 77 - - - 8 - Table 4.5 Tasks parameters in the four scenarios. Task (coordinate x, coordinate y)ϕjψj(s) S-I S-II S-III S-IV S-I S-II S-III S-IV S-I S-II S-III S-IV t1(392.5, 302) (227.5, 176) (821.5, 218) (161.5, 218) 3 2 5 8 10 40 0 30 t2(458.5, 302) (689.5, 472) (1085.5, 134) (29.5, 176) 1 3 2 6 10 0 0 23 t3(1151.5, 472) (491.5, 134) (1019.5, 48) (524.5, 344) 5 5 2 3 30 20 40 69 t4(623.5, 176) (623.5, 5) (1151.5, 48) (788.5, 5) 1 4 4 1 20 40 60 48 t5(788.5, 344) (1085.5, 344) (986.5, 218) (161.5, 91) 2 2 5 7 20 70 20 93 t6(854.5, 48) (359.5, 48) (920.5, 5) (425.5, 134) 1 2 4 2 10 80 25 109 t7- (62.5, 429) (1052.5, 91) (491.5, 91) - 10 2 8 - 30 30 57 t8- (656.5, 472) (1151.5, 5) (788.5, 429) - 2 1 6 - 50 10 21 t9- - (1118.5, 218) (161.5, 302) - - 4 1 - - 20 107 t10 - - - (887.5, 302) - - - 6 - - - 68 t11 - - - (326.5, 48) - - - 3 - - - 21 t12 - - - (590.5, 5) - - - 3 - - - 23 t13 - - - (920.5, 344) - - - 6 - - - 38 t14 - - - (1052.5, 386) - - - 3 - - - 11 t15 - - - (260.5, 48) - - - 1 - - - 109 Table 4.6 GA tuning parameters. PS FEFCFM1 FM2 Max Gen. Max Stall Gen. 100 0.1 0.7 0.5 0.5 50 30 4.7 Results This section shows and discusses the simulation results of the application of the proposed algorithms for the proposed scenarios. 4.7.1 Scenario I In this case, the best allocation obtained by solving Scenario I using the GA, which is also the most frequently obtained allocation among the of the GA, is the optimal allocation, obtained with the BaB and can be seen in Fig. 4.11 (where ground robots appear in purple while the aerial robots appear in blue). The optimal allocation implies that r1 and r2 stay still, r3 performs t1 and t2 sequentially, and r4 performs the rest of the tasks in the following order: t3,t5,t6and t4. The best initial solution obtained for Scenario I, which corresponds to the assignment problem for distances that was depicted in eq. (4.16) , is shown in Fig. 4.12. As it can 5 Game Theory Based Clustering for MRTA Bearing in mind the application mentioned in Chapter 2, we consider here the problem of assigning a set of measuring tasks, presumably generated by an upper layer of the control system, to the set of available robots. As we saw in Chapter 4 the combinatorial nature of the problem can render its computation infeasible from a computational viewpoint, especially for large-scale systems where many robots and tasks can be expected. To relieve this issue, we propose a clustering approach to group sets of robots and tasks, which will now become the players of a game. Clustering is a halfway approach between centralized and distributed methods, which can relieve computation burden. In this way, sets of loosely coupled tasks and robots are divided into clusters, which can be dealt with in parallel. Indeed, clustering has been applied to the formation of robot coalitions to perform complex tasks, e.g., in (Shehory and Kraus, 1998a). Also, in (Dutta et al., 2019), a correlation clustering technique that enables similar robots to form coalitions is explored; and, in (Arif, 2021), an evolutionary algorithm is proposed for the case in which robots must participate in different clusters. However, splitting big MRTA problems into several smaller, manageable and mutually exclusive problems is an NP-hard problems, e.g., the capacitated clustering problem (Mulvey and Beck, 1984; Negreiros and Palhano, 2006), which consists of dividing a set of elements into clusters with limited capacity and maximum similarity within the clusters, is illustrative of the difficulties in this regard. Therefore, heuristic methods such as (Brimberg et al., 2019) become necessary. Clustering fits naturally in the cooperative game theory framework, where intelligent rational decision-makers, known as players, are involved in cooperative situations or games. In this context, players can make sacrifices in terms of their own welfare to improve overall performance, and any group of players or cluster is free to make agreements and unify strategies (Peleg and Sudhölter, 2007). Note that the application of cooperative game theory to engineering problems is common, with contributions in the field of cooperative structures (Gilles, 2010), communication networks (Chasparis and Shamma, 2013; Han et al., 2012; Saad et al., 2009; Zhang et al., 2011), smart grids (Marzband et al., 2017; 65 66 Chapter 5. Game Theory Based Clustering for MRTA Saad et al., 2012), MRS (Jaleel and Shamma, 2020), and also MRTA problems (Arslan et al., 2007; Chapman et al., 2009). The proposed method is based on the ordering of players introduced by the Shapley value (Shapley, 1953), the best-known solution concept in cooperative game theory, with applications in water systems (Muros et al., 2018), social networks (Narayanam and Narahari, 2011), wine ranking (Ginsburgh and Zang, 2012a), and pollution reduction (Petrosjan and Zaccour, 2003), among others. Furthermore, the Shapley value has been used in multiagent systems to form coalitions due to its relationship with the marginal contribution of each agent (Liao et al., 2014). Also, the Shapley value has been widely used in the literature to perform rankings, e.g., (Ginsburgh and Zang, 2012b; Maestre et al., 2014; Muros et al., 2017b; Narayanam and Narahari, 2011; Zhang et al., 2014). Finally, a key feature in this context is that the Shapley value can be approximated in polynomial time using randomized methods (Castro et al., 2009; Tarkowski et al., 2018) so as to relieve the computation burden of the Shapley value. In particular, robots and tasks are considered as players of a game and are ranked according to their Shapley value, which captures their relevance in the MRTA scenario. Based on this value, an algorithm is proposed to group the different players into clusters. In this way, and as it will be shown in the results section, our method can outperform other heuristics such as genetic algorithms, even for mild timing constraints regarding the computation of the MRTA problem. Moreover, due to the relevance of the Shapley value in machine learning applications, where it is used to assess the quality of data points (Fryer et al., 2021; Ghorbani and Zou, 2019; Pradhan et al., 2022), new techniques are being developed that can deal with thousands of data points/players (Ghorbani et al., 2020; Jia et al., 2019; Kwon et al., 2021). The rest of this Chapter is organized as follows. The problem statement is introduced in Section 5.1. Some game theory concepts used in this work are presented in Section 5.2. In Section 5.3, the Shapley value-based clustering algorithm for MRTA problems is proposed. Likewise, in Section 5.4, a large-scale case study is introduced to assess the algorithm, which is also compared to a GA and a fully centralized approach. Finally, conclusions are given in Section 5.5. 5.1 Problem Statement We consider the sets R and T mentioned in Chapter 2.1. For convenience, let us define set P=R∪T that groups robots and tasks. To solve the MRTA problem, the elements inside R and T need to be properly matched considering that, as in Chapter 4, we focus on ST-SR-TA MRTA problems (Gerkey and Matarić, 2004; Gerkey and Mataric, 2004; Korsah et al., 2013). Also, as in previous chapters and for the sake of simplicity, robots are governed by a controller that reaches the assigned task positions, avoiding static and dynamic obstacles (Čáp et al., 2015; Zhu, 1991). Let us introduce the MRTA clustering as G={G1,G2,...,GNc},(5.1) where Gh,h=1,...,Ncare disjoint clusters of robots and tasks verifying 5.1 Problem Statement 67 [ h Gh=P,Gha∩Ghb=/0,∀Gha,Ghb∈G,Gha=Ghb.(5.2) Note that the number of different clustering alternatives for G is upper bounded by the Bell Number (Becker and Riordan, 1948), BP= P ∑ s=0 1 s! s ∑ j=0 (−1)s−js jjP!,(5.3) with P=|P| . In particular, clusters Gh composed only by robots or tasks are allowed in this definition. Nevertheless, as it will be shown later, they will be discouraged via penalties, specially those containing only tasks because task-only clusters lead to tasks not being performed because tasks can only be executed by robots. Hence, in clusters without robots, tasks would remain undone. 5.1.1 MRTA Objective Function Centralized approaches for MRTA problems rely on functions that evaluate allocations, i.e., the plan to perform the tasks, according to criteria such as the distance traveled by robots, their battery level, and the time to accomplish the missions (Jones et al., 2011a; Jose and Pratihar, 2016; Martin et al., 2021a; Panchu K. et al., 2018; Tolmidis and Petrou, 2013). Following the work developed in (Martin et al., 2021a) and depicted in Chapter 4, the multi-criteria evaluation function used is (4.3) , which includes the considerations listed in Chapter 2.7. 5.1.2 Cluster Evaluation Since it may not be possible to solve (4.3) in real time for large problems, we propose to divide robots and tasks into clusters or coalitions, where the latter is the term employed in game theory to refer to subset of players. To this end, we introduce a function JGP(GP) , to evaluate a coalition GP⊆Paccording to its robots and tasks, considering: • The average distance between robots and tasks so that larger distances are penalized. •The velocity of the robots in the cluster. The faster the robots, the better. •The battery level of the robots in the cluster. Again, the higher it is, the better. • The penalties of robots and tasks in (4.3) should be accounted for, so that neither critical tasks nor penalized robots are concentrated in a single cluster. • The operation time to accomplish since τj(γ) in (4.3) depends on these operation times and also on the travel times spent by the robots. •Large cardinality of GP,|GP|=G, should also be penalized. According to all these premises, we define 68 Chapter 5. Game Theory Based Clustering for MRTA JGP(GP) = ζ1(GP) ∑ i⊆GP JGP i−α1·¯ D(GP)−α2·G+¯ JGP!,(5.4) with JGP i=(α3·bi+α4·vi−α5·λi,if i∈R, −(α6·ϕi+α7·ψi),if i∈T,(5.5) and ζ1(GP) = (0,if {GP⊆R}∨{GP⊆T}∨{GP≡/0}, 1,otherwise.(5.6) Several terms are considered by cost function (5.4) . First, JGP i is related to the cost of each element i∈P , while ¯ D(GP) is the average distance from robots to tasks in coalition GP . Finally, term ζ1(GP) filters the cases where there is only one type of player inside coalition GP , i.e., either robots or tasks, with ζ1(GP) = 0 in those coalitions, and ¯ JGP is an arbitrarily large number that guarantees that JGP(GP) is positive. Recall λi and ϕi were previously defined in (4.3) . Note that α1 to α7 are nonnegative weighting scalars with different effects over the resulting clustering. Using larger α1 reduces the admissible distance between robots and tasks; similarly, increasing α2 decreases the number of players in the formed clusters; also, increasing α3 , α4 , α5 , α6 , and α7 prioritizes aspects such as the battery level, velocity, robot penalty, task relevance, and operation time, respectively. These weights can be set to 1 initially and then modified on the basis of the system designer’s experience and their impact should be simulated in a trial and error fashion. As stated, ¯ JGP balances the computation to ensure JGP(GP)≥0 . To this end, ¯ JGP is pre-computed as ¯ JGP=ˆ JGP(GP w), where ˆ JGP(GP) = ζ1(GP) ∑ i⊆GP JGP i−α1·¯ D(GP)−α2·G!.(5.7) Note that (5.7) is equivalent to (5.4) . Even when it is not possible to know beforehand the worst scenario it is possible to find clear candidates for GP w taking into account that, in general, robots add value to the output, and tasks subtract it: • In the set of coalitions with one robot and all tasks, i.e., the R coalitions that verify GP=T∪riwith ri∈R. • For larger α1 , the worst coalition GP w might be composed of one robot and one task, i.e., the T·R coalitions that verify G=2 , GP⊈R , GP⊈T . Note that when there are only one robot and one task very far away one from the other, ¯ D(GP) becomes very relevant for the function output. •In case α2is dominant in (5.4), the set Pmight become the worst coalition GP w. 5.1 Problem Statement 69 Therefore, by pre-computing the aforementioned R·(1+T) + 1 possibilities, it is possible to find a suitable candidate for ¯ JGP . Nevertheless, note that due to parameter α5 , robots do not always add value to the output of (5.4) , so it is not possible to assure that this initial ¯ JGP corresponds to the worst possible coalition. For these cases, ¯ JGP needs to be updated if a negative value of JGP(GP)is detected. Note that: •In contrast to other works, we consider that tasks, which are characterized by their location and other parameters such as their relevance and and the operation time required to perform them, can also be players in the cooperative game and also have an effect on the assignment cost function (4.3). • The MRTA problem (4.3) has inspired the design of function (5.4) , but alternative evaluation functions such as the one presented in (Zhang and Parker, 2013), and different MRTA objective functions could be used as well. • Parameter α2 in (5.4) penalizes the coalition size. Thus, lower/higher values of α2 will favor larger/smaller coalitions. Note that the cluster size affects the allocation computation by (4.3) . Indeed, larger clusters, gain optimality at the expense of computation cost. 5.1.3 Clustering Evaluation We define a function V(G)to assess clustering Gas V(G) = ζ2(G)∑ Gh⊆G JGP(Gh), s.t. (5.2), (5.8) with JGP(·) given by (5.4) , and where ζ2(G) = 0 if any Gh∈G is composed only by tasks and is set to 1 otherwise. This function V(G) can be maximized to obtain the optimal allocation. In the following section, a clustering method based on game theory is proposed, where it will be necessary to calculate the value of all possible coalitions (5.4) by means of GP⊆P . Then, the cost for each specific clustering solution G(5.1) is evaluated by (5.8) . Note that the cost of clusters Gh∈G is computed by function J in (5.8) since these clusters belong to the full set of coalitions, i.e., Gh⊆P. For instance, let us consider a set P composed by R=5 robots and T=12 tasks, i.e., P=17 players, as shown in Fig. 5.1 and Fig. 5.2. The cost of the corresponding 217 coalitions can be computed by (5.4) . Then, a possible clustering G={G1,G2,G3},Gh⊆P , the cost of which would be computed by (5.8) , and the parameters taken into account are also shown in Fig. 5.3. 5.1.4 Proposed Control Scheme To conclude this section, the overall control scheme proposed in this work is illustrated in Fig. 5.4 in which the number of elements in each cluster is respectively denoted by N1,N2,...,NNc . On the top, a master control layer is in charge of performing the clustering, 70 Chapter 5. Game Theory Based Clustering for MRTA 500 m 500 m 1 1 7 8 2 5 3 4 6 10 9 11 12 3 2 4 5 6 16 10 12 9 8 1 3 7 2 15 14 13 11 17 4 5 t t t t t t t t t t t t r r r r r Figure 5.1 17-agent system example. 500 m 500 m φ12 ψ12 φ11 ψ11 φ8 ψ8 φ9 ψ9 φ7 ψ7 φ10 ψ10 φ6 ψ6 φ5 ψ5 φ4 ψ4 φ3 ψ3 φ2 ψ2 φ1 ψ1 b1 v1 λ1 b3 v3 λ3 b4 v4 λ4 b5 v5 λ5 b2 v2 λ2 Figure 5.2 Parameters of robots and tasks in (4.3). that is, dividing the robots and tasks into Nc different groups (recall (5.1) ), to minimize function (5.8) . Once the clustering is performed, a first intermediate layer solves the MRTA problem for each cluster Gh,h=1,...,Nc , providing a suitable allocation by establishing the sequence of tasks performed by each robot. Note that this allocation problem for 5.2 Game Theory Viewpoint 71 500 m 500 m φ12 ψ1 φ11 ψ1 φ8 ψ8 φ9 ψ9 φ7 ψ7 φ10 ψ1 0 φ6 ψ φ5 ψ5 φ4 ψ4 φ3 ψ3 φ2 ψ2 φ1 ψ1 b1 v1 λ1 b3 v3 λ3 b4 v4 λ4 b5 v5 λ5 b2 v2 λ2 1 1 7 8 2 5 3 4 6 10 9 11 12 3 2 5 6 16 10 12 9 8 1 3 7 2 15 14 13 11 17 5 t t t t t t t t t t t t r r 4 4 r r r Figure 5.3 Clustering based on parameters in (4.3). each cluster could be solved either in a centralized (Al-Yafi et al., 2009; Brumitt and Stentz, 1998; Coltin and Veloso, 2010) or a distributed (Choi et al., 2009; Gao et al., 2009; Giordani et al., 2010) fashion. In this chapter, we consider the centralized solution introduced in (Martin et al., 2021a) that was depicted in Chapter 4, which is based on solving (4.3) . Finally, once the allocation is sent to the robots, two individual bottom layers are also locally required. The high-level layer plans the trajectory of each robot to perform its tasks and avoid obstacles. Also, a low-level layer manages technical aspects such as battery level, data buffer capacity, etc. 5.2 Game Theory Viewpoint In this work, the pair (P,J J JGP) will be interpreted as a cooperative game with transferable utility (TU-game), where P is the set of players and the cost function J J JGP assigns cost JGP(GP) defined by (5.4) to each coalition GP⊆P . Note that a necessary condition for TU-games is JGP(/0) = 0 , which is assured by ζ1(GP) in (5.4). From the different cooperative game theory payoff rules available in the literature, we will consider the Shapley value (Shapley, 1953) to compute the relevance of the coalitions among the different players – robots and tasks – involved in our game. The Shapley value assigns vector φ φ φ(P,J J JGP) = hφi(P,J J JGP)ii∈Pto game (P,J J JGP), with φi(P,J J JGP) = ∑ GP⊆P:i/∈GP G!(P−G−1)! P! [JGP(GP∪{i})−JGP(GP)].(5.9) 72 Chapter 5. Game Theory Based Clustering for MRTA MASTER CONTROL LAYER (PERFORMS THE CLUSTERING) MRTA FOR CLUSTER 1 ... ROBOT 1 HIGH-LEVEL CONTROL MRTA FOR CLUSTER 2 MRTA FOR CLUSTER NC ... ROBOT N1 HIGH-LEVEL CONTROL ... ROBOT 1 HIGH-LEVEL CONTROL ROBOT N2 HIGH-LEVEL CONTROL ... ROBOT 1 HIGH-LEVEL CONTROL ROBOT NNc HIGH-LEVEL CONTROL ... ROBOT 1 LOW-LEVEL CONTROL ROBOT N1 LOW-LEVEL CONTROL ... ROBOT 1 LOW-LEVEL CONTROL ROBOT N2 LOW-LEVEL CONTROL ... ROBOT 1 LOW-LEVEL CONTROL ROBOT NNc LOW-LEVEL CONTROL LOCAL CENTR. DISTR./CENTR. Figure 5.4 Control scheme of the proposed approach. Here, J(GP∪{i})−JGP(GP) = µi(P,J J JGP) is the marginal contribution of player g when it incorporates to coalition GP and G!(P−G−1)! P! is the probability for player g to join GP . That is, the Shapley value can be interpreted as the expected marginal contribution for each player when joining a coalition randomly. The key idea of our proposal is to use the ordering of players generated by the Shapley value to perform the clustering, as it has been previously done in other works for coalition formation (Liao et al., 2014; Shehory and Kraus, 1998b). Even when this value is oriented to the distribution of payoffs among players in the grand coalition, the corresponding payoff is based on their average marginal contribution when they join a random coalition. Therefore, the Shapley value can be interpreted as a measure of the relevance of every player that is weighted across the set of possible coalitions where it can participate, which is why it can be used as a ranking criterion for the players. Note that the Shapley value was introduced axiomatically as the only solution concept that satisfies the properties of null player,symmetry,additivity and efficiency (Shapley, 1953). In particular, the efficiency property states that the sum of all Shapley values remains constant by satisfying ∑ i∈P φi(P,J J JGP) = J(P),(5.10) which will be of particular interest in the clustering algorithm presented below. Alternatively to (5.9), the Shapley value can be rewritten in terms of all possible 5.2 Game Theory Viewpoint 73 orderings of players in P coming into a coalition, i.e., P! . Hence, assuming equiprobable orderings, the Shapley value can be computed by (Weber, 1988) φi(P,J J JGP) = 1 P! ∑ π∈Π(P) µπ i(P,J J JGP),∀i∈P,(5.11) , where µπ i(P,J J JGP) = J({j∈P|π(j)≤π(i)}) −J({j∈P|π(j)<π(i)}),(5.12) is the marginal contribution of player g to the players that are ranked before it in permutation π, and with Π(P)being the collection of all permutations. Following (Xu et al., 2008), it is possible to find a matrix expression for the Shapley value. Consider a matrix L∈RP×2P , where the rows correspond to each player g∈P and the columns to the different coalitions GP⊆P in the lexicographic order. Given a game (P,J J JGP) , the Shapley value φ φ φ(P,J J JGP) can be represented by Shapley standard matrix L as φ φ φ(P,J J JGP) =      φ1 φ2 . . . φP     =L·         JGP(GP 0) JGP(GP 1) JGP(GP 2) . . . J(GP 2P−1)         =L·J J JGP,(5.13) with each component of J J JGP given by (5.4), and where each element of L , denoted by liGP , is defined as liGP=   (G−1)!(P−G)! P! i∈GP, −G!(P−G−1)! P! i/∈GP. (5.14) Note that it is possible to rewrite the Shapley value of every player i∈GP in terms of (5.14) by the following expression: φi(P,J J JGP) = ∑ GP⊆P liGP·JGP(GP).(5.15) Coalitions composed solely of either robots or targets are excluded by means of term ζ1(GP) . Nevertheless, it is possible to use other methods in the line of (Hiller, 2016, 2018) to exclude prohibited coalitions via a Shapley value redefinition. 5.2.1 Estimating the Shapley Value Computing the Shapley value by (5.9) becomes increasingly difficult due to the exponential growth of the problem size as it can be seen in Fig. 5.5. 80 Chapter 5. Game Theory Based Clustering for MRTA 0 50 100 150 200 250 300 350 400 450 500 x coordinate (m) 0 50 100 150 200 250 300 350 400 450 500 y coordinate (m) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 1 2 3 45 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 Robot Task Figure 5.7 Randomly generated scenario (P=60). Table 5.4 Randomly generated scenario parameters for robots and tasks. ribiviλitjψjϕjtjψjϕjtjψjϕj r197 6 5 t18 3 t16 5 1 t31 4 2 r299 8 1 t21 5 t17 9 5 t32 8 2 r343 18 1 t38 4 t18 5 1 t33 10 4 r490 11 4 t43 3 t19 6 4 t34 4 3 r52 8 3 t59 2 t20 5 1 t35 3 4 r61 18 1 t68 4 t21 6 2 t36 7 5 r76 20 2 t75 3 t22 9 4 t37 10 3 r858 8 4 t86 3 t23 2 2 t38 3 3 r98 13 5 t96 3 t24 9 2 t39 5 1 r10 80 8 3 t10 2 2 t25 1 3 t40 2 1 r11 70 5 5 t11 8 4 t26 1 5 t41 10 1 r12 86 8 1 t12 5 3 t27 6 5 t42 8 5 r13 31 10 4 t13 3 5 t28 3 1 t43 6 4 r14 58 6 5 t14 7 3 t29 8 2 t44 8 3 r15 62 5 3 t15 9 4 t30 2 4 t45 6 1 clustering results of Fig. 5.9, where various colors represent the 10 different clusters while diverse line styles refer to different robot trajectories in the same cluster. 5.4 Case Study 81 0 50 100 150 200 250 300 350 400 450 500 x coordinate (m) 0 50 100 150 200 250 300 350 400 450 500 y coordinate (m) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 1 2 3 45 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 Robot Task Figure 5.8 Considered scenario solved via a centralized allocation. 0 50 100 150 200 250 300 350 400 450 500 x coordinate (m) 0 50 100 150 200 250 300 350 400 450 500 y coordinate (m) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 1 2 3 45 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 Figure 5.9 Scenario clustered by the SVCA. Since Fig. 5.9 its too messy, we have split the different clusters in different figures: In Fig. 5.10 cluster G1 both the allocation obtained in this cluster are shown, in Fig. 5.11 is equivalent for cluster G2 and the allocation obtained in this cluster are shown, in Fig. 5.12 is equivalent for cluster G3 , in Fig. 5.13 is equivalent for cluster G4 , in Fig. 5.14 is equivalent for cluster G5 , in Fig. 5.15 is equivalent for cluster G6 , in Fig. 5.16 is equivalent for cluster G7 , in Fig. 5.17 is equivalent for cluster G8 , in Fig. 5.18 cluster G9 , in Fig. 5.19 is equivalent for cluster G10 , in Fig. 5.20 is equivalent for cluster G11 , in Fig. 5.21 is equivalent for cluster G12, and finally, Fig. 5.22 is equivalent for cluster G13. 82 Chapter 5. Game Theory Based Clustering for MRTA 0 50 100 150 200 250 300 350 400 450 500 x coordinate (m) 0 50 100 150 200 250 300 350 400 450 500 y coordinate (m) 1 11 15 20 21 27 (a) Cluster G1. 0 50 100 150 200 250 300 350 400 450 500 x coordinate (m) 0 50 100 150 200 250 300 350 400 450 500 y coordinate (m) 1 11 15 20 21 27 (b) Cluster G1solved. Figure 5.10 Cluster G1in SVCA. 0 50 100 150 200 250 300 350 400 450 500 x coordinate (m) 0 50 100 150 200 250 300 350 400 450 500 y coordinate (m) 2 2 6 32 37 45 (a) Cluster G2. 0 50 100 150 200 250 300 350 400 450 500 x coordinate (m) 0 50 100 150 200 250 300 350 400 450 500 y coordinate (m) 2 2 6 32 37 45 (b) Cluster G2solved. Figure 5.11 Cluster G2in SVCA. 5.4 Case Study 83 0 50 100 150 200 250 300 350 400 450 500 x coordinate (m) 0 50 100 150 200 250 300 350 400 450 500 y coordinate (m) 3 36 39 (a) Cluster G3. 0 50 100 150 200 250 300 350 400 450 500 x coordinate (m) 0 50 100 150 200 250 300 350 400 450 500 y coordinate (m) 3 36 39 (b) Cluster G3solved. Figure 5.12 Cluster G3in SVCA. 0 50 100 150 200 250 300 350 400 450 500 x coordinate (m) 0 50 100 150 200 250 300 350 400 450 500 y coordinate (m) 4 3 4 9 38 41 (a) Cluster G4. 0 50 100 150 200 250 300 350 400 450 500 x coordinate (m) 0 50 100 150 200 250 300 350 400 450 500 y coordinate (m) 4 3 4 9 38 41 (b) Cluster G4solved. Figure 5.13 Cluster G4in SVCA. 84 Chapter 5. Game Theory Based Clustering for MRTA 0 50 100 150 200 250 300 350 400 450 500 x coordinate (m) 0 50 100 150 200 250 300 350 400 450 500 y coordinate (m) 5 6 22 33 40 44 (a) Cluster G5. 0 50 100 150 200 250 300 350 400 450 500 x coordinate (m) 0 50 100 150 200 250 300 350 400 450 500 y coordinate (m) 5 6 22 33 40 44 (b) Cluster G5solved. Figure 5.14 Cluster G5in SVCA. 0 50 100 150 200 250 300 350 400 450 500 x coordinate (m) 0 50 100 150 200 250 300 350 400 450 500 y coordinate (m) 7 9 14 18 (a) Cluster G6. 0 50 100 150 200 250 300 350 400 450 500 x coordinate (m) 0 50 100 150 200 250 300 350 400 450 500 y coordinate (m) 7 9 14 18 (b) Cluster G6solved. Figure 5.15 Cluster G6in SVCA. 5.4 Case Study 85 0 50 100 150 200 250 300 350 400 450 500 x coordinate (m) 0 50 100 150 200 250 300 350 400 450 500 y coordinate (m) 8 1 8 29 (a) Cluster G7. 0 50 100 150 200 250 300 350 400 450 500 x coordinate (m) 0 50 100 150 200 250 300 350 400 450 500 y coordinate (m) 8 1 8 29 (b) Cluster G7solved. Figure 5.16 Cluster G7in SVCA. 0 50 100 150 200 250 300 350 400 450 500 x coordinate (m) 0 50 100 150 200 250 300 350 400 450 500 y coordinate (m) 10 19 28 34 35 (a) Cluster G8. 0 50 100 150 200 250 300 350 400 450 500 x coordinate (m) 0 50 100 150 200 250 300 350 400 450 500 y coordinate (m) 10 19 28 34 35 (b) Cluster G8solved. Figure 5.17 Cluster G8in SVCA. 86 Chapter 5. Game Theory Based Clustering for MRTA 0 50 100 150 200 250 300 350 400 450 500 x coordinate (m) 0 50 100 150 200 250 300 350 400 450 500 y coordinate (m) 11 5 7 26 (a) Cluster G9. 0 50 100 150 200 250 300 350 400 450 500 x coordinate (m) 0 50 100 150 200 250 300 350 400 450 500 y coordinate (m) 11 5 7 26 (b) Cluster G9solved. Figure 5.18 Cluster G9in SVCA. 0 50 100 150 200 250 300 350 400 450 500 x coordinate (m) 0 50 100 150 200 250 300 350 400 450 500 y coordinate (m) 12 10 12 42 43 (a) Cluster G10. 0 50 100 150 200 250 300 350 400 450 500 x coordinate (m) 0 50 100 150 200 250 300 350 400 450 500 y coordinate (m) 12 10 12 42 43 (b) Cluster G10 solved. Figure 5.19 Cluster G10 in SVCA. 5.4 Case Study 87 0 50 100 150 200 250 300 350 400 450 500 x coordinate (m) 0 50 100 150 200 250 300 350 400 450 500 y coordinate (m) 13 24 31 (a) Cluster G11. 0 50 100 150 200 250 300 350 400 450 500 x coordinate (m) 0 50 100 150 200 250 300 350 400 450 500 y coordinate (m) 13 24 31 (b) Cluster G11 solved. Figure 5.20 Cluster G11 in SVCA. 0 50 100 150 200 250 300 350 400 450 500 x coordinate (m) 0 50 100 150 200 250 300 350 400 450 500 y coordinate (m) 14 13 17 25 (a) Cluster G12. 0 50 100 150 200 250 300 350 400 450 500 x coordinate (m) 0 50 100 150 200 250 300 350 400 450 500 y coordinate (m) 14 13 17 25 (b) Cluster G12 solved. Figure 5.21 Cluster G12 in SVCA. 88 Chapter 5. Game Theory Based Clustering for MRTA 0 50 100 150 200 250 300 350 400 450 500 x coordinate (m) 0 50 100 150 200 250 300 350 400 450 500 y coordinate (m) 15 16 23 30 (a) Cluster G13. 0 50 100 150 200 250 300 350 400 450 500 x coordinate (m) 0 50 100 150 200 250 300 350 400 450 500 y coordinate (m) 15 16 23 30 (b) Cluster G13 solved. Figure 5.22 Cluster G13 in SVCA. The specific distance traveled by the robots and the tasks completion time due to the clustering solution is shown in Table 5.6 and the clusters formed can be consulted in Table 5.7. The clustering performance is given by V(G) = 10921.20 , while the allocation cost computed by (4.3) is J=83392.52 . Likewise, the overall time to perform the clustering and achieve the allocation is tc clu +tc all =8.53+1.57 =10.10 s. Note that the computation time is significantly reduced with respect to the centralized solution, though at the cost of a performance loss. 5.4.3 Genetic Algorithms (GA)-Based Clustering Since problem (5.8) is highly nonlinear, it is not possible to find the optimal clustering without evaluating all the possibilities, which may not be computationally feasible. For this reason, in order to measure the performance of the proposed algorithm, we have considered the well-known Genetic Algorithm (GA) (Kudova, 2007; Maulik and Bandyopadhyay, 2000; Mitchell, 1998). The rationale of GA is to start from a random population of genes representing possible solutions and improve them by mixing and mutating them. The gen used in GA has been defined as illustrated in Table 5.5, and it can be expressed mathematically as a vector g g g=ggg∈P , where gg∈[1,N] is the cluster containing player g . Note that the number of genes generated increases fast because the number of variations with repetitions is given by PP . Finally, it is interesting to introduce the population size, i.e., the amount of genes that the algorithm uses in each generation. Indeed, increasing this parameter can improve the optimality of the solution achieved, but it also increases the overall computation time. Table 5.5 Genome involved in the GA. Player 1··· g··· P Cluster [1,P]··· [1,P]··· [1,P] 5.4 Case Study 89 In this work, the GA has been implemented using ga in Matlab ® with default parameters. Specifically, the same scenario with P=60 as before has been solved using GA. The clustering results are shown in Fig. 5.23, where, as in the previous case, different colors represent the eight resulting clusters and various line styles refer to different robot trajectories in the same cluster. 0 50 100 150 200 250 300 350 400 450 500 x coordinate (m) 0 50 100 150 200 250 300 350 400 450 500 y coordinate (m) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 Figure 5.23 Scenario clustered by the GA. Similarly to what happened in SVCA, since Fig. 5.23 is messy, we have splitted the figure into several figures where Fig. 5.24 correspond to cluster G1 , Fig. 5.25 to cluster G2 , Fig. 5.26 to cluster G3 , Fig. 5.27 to cluster G4 , Fig. 5.28 to cluster G5 , Fig. 5.29 to cluster G6, Fig. 5.30 to cluster G7, and Fig. 5.31 to cluster G8. 0 50 100 150 200 250 300 350 400 450 500 x coordinate (m) 0 50 100 150 200 250 300 350 400 450 500 y coordinate (m) 1 8 1 5 12 14 29 40 42 (a) Cluster G1. 0 50 100 150 200 250 300 350 400 450 500 x coordinate (m) 0 50 100 150 200 250 300 350 400 450 500 y coordinate (m) 1 8 1 5 12 14 29 40 42 (b) Cluster G1solved. Figure 5.24 Cluster G1in GA. 96 Chapter 5. Game Theory Based Clustering for MRTA 0 5 10 15 20 25 30 35 3000 3500 4000 4500 5000 5500 6000 6500 7000 7500 SVCA GA Figure 5.33 V(G)vs. Computation time (s). and achieving a feasible solution much faster than centralized schemes. It is important to remark that the proposed method does not only provide a clustering, but also qualitative information regarding how useful a robot is or how demanding a task is. In the current implementation, the proposed method can manage problems of about hundreds of agents. Therefore, there might be limitations in the number of players derived from the available time for the computations. Nevertheless, parallel computing techniques can be implemented to reduce the calculation time because the algorithm in (Castro et al., 2009, 2017) allows to compute each sample independently. Furthermore, the new techniques being developed to compute the Shapley value allow dealing with problems with thousands of players, e.g., in machine learning applications and might enhance the applicability of our method. 6 Predictive MRTA with dynamic tasks In previous chapters tasks were considered static. In this chapter we will consider a different approach with tasks moving across the solar plant since TA MRTA problems become more challenging when tasks move (Jin and Li, 2016; Jin et al., 2019). An algorithm that predicts the evolution of the tasks and robots as it is done in the model predictive control (MPC) framework is proposed. Predictive control methods in the MRS field have been proposed in other works, e.g., (Kuhn et al., 2011), where a decentralized model-predictive control strategy is used to address the observation of multiple moving targets; (Augusto de Holanda et al., 2020), where a non-linear model predictive control strategy is used to dynamically set the formation leader, and (Gaggero et al., 2019), where it is used to compute the time to change the assignment of tasks to robots. The contents of this chapter have been first published in (Martin et al., 2022), and further developed in (Martin et al., 2023b). The proposed algorithm predicts the evolution of tasks and robots performs calculation over an allocation horizon and applies the first allocated task to each robot. Unlike other algorithms based on solving a non-convex optimization problem, this method is based on transforming the problem into a Linear Program (LP) one, which is known to have a much lower computation burden. The rest of this chapter is organized as follows: Section 6.1 presents the problem statement and assumptions; in Section 6.2, the proposed algorithm for the spatial estimation is depicted; in Section 6.3, the simulation environment and results are shown and discussed. Finally, conclusions and future investigation lines are given in Section 6.4. 6.1 Problem Statement In this work qj∈R2 corresponds to the position of task tj , and not to the central position of a certain cell, cj . This way, the velocity of task tj is denoted as ˙qj∈R2 , with |˙qj|>0 . Also, for the sake of simplicity, let us consider OG=/0 and OA=/0 . An example can be seen in Fig. 6.1, where distances from robots to task t1 are represented in red, to task t2 in 97 98 Chapter 6. Predictive MRTA with dynamic tasks blue, and to t3in pink. Also, task trajectories can be seen in green. 01234567 X-axis (m) -50 0 50 100 150 200 Y-axis (m) Figure 6.1 Example with 3robots and 3tasks. Let us define set K={1,2,...,K}⊂N as the allocation horizon, being K its length. Then, let δi jk ∈{0,1} take value 1 when robot ri is allocated to fulfill task tj as its k -th mission and 0 otherwise, with k∈K . Note that there is not a straight forward relation between kand the time when a robot will perform its allocated tasks. An example of the described variables can be seen in Fig. 6.2, where tasks are considered static. Note that blue circles represent robots and red circles represent tasks. Robot r1 performs first task t1 and then task t2 ; robot t2 performs task t3 . In this example, task t1 completion time, τ1 , is the time it takes for robot r1 to reach t1 plus ψ1 ; task t2 completion time, τ2 , is τ1 , plus the time it takes robot r1 to go from task t1 to task t2 and perform task t2 ; and task t3 completion time, τ3 , is the time it takes for robot r2 to reach task t3 plus ψ3 . Likewise, the distance traveled by r1 is the distance from its initial position to t1 plus the distance from t1to t2; and the distance traveled by r2is that from its initial position to t3. Finally, let us also define δi= [δi11δi12 ···δi1Kδi21 ···δiTK] as a vector representing the possible assignments of robot ri and δ= [δ1···δR] as the vector representing the assignment of all the robots. Let us consider the cost function J described in Chapter 4.2. When computing an allocation, the positions of robots, tasks, and obstacles in the plant are known. Therefore, the distance that robot ri has to travel to reach task tj , i.e., di j (as shown in Fig. 6.1 for a case with no obstacles) can be computed. Likewise, the time that the robot rineeds to reach the task tj,ti j, can also be computed as τi j =di j vi .(6.1) Furthermore, since the velocity of robots and tasks is known, we can compute an 6.1 Problem Statement 99 0.5 1 1.5 2 2.5 X-axis (m) 0.5 1 1.5 2 2.5 3 3.5 Y-axis (m) Figure 6.2 Example of variables. 0 2 4 6 X-axis (m) 1 2 3 4 5 6 Y-axis (m) Figure 6.3 Interception point between robot riand task tj,ℓi j. interception point between them, ℓi j ∈R2 , to improve the computed distances and times as it is illustrated in Fig. 6.3 for the case of |¨qj|=0 and with no obstacles. τℓ i j is the time that it takes for robot ri to intercept task tj . In case that task tj is allocated to robot ri , the robot will move towards ℓi j at speed viconsidering the movement of the task. Taking into account that qj(t) and qqj(t=0) are known, we can estimate the location of task tjin kas 100 Chapter 6. Predictive MRTA with dynamic tasks q′jk =qj(τA (k−1)).(6.2) with q′jk being the estimated position of task tj when starting its allocated k -th task, and τA (k−1) the elapsed time from the moment the allocation is performed to the moment when robots start their k-th task. However, predicting the location of the robots is challenging because it depends highly on the allocation. Moreover, robots move towards their previously allocated active task. In (Martin et al., 2022), this issue was simplified by predicting the movement of tasks and considering that robots move towards their gravity center. However, this estimation differs rapidly from reality. To improve it, in (Martin et al., 2023b) we proposed to estimate the position of tasks and compute the distances between them, since we assume that in the next step k+1 the robot will have moved to its interception point with one of its tasks. This way we have dl i jk as the distance from robot ri to its interception point with task tj as its k-th mission if task tlwas done previously. di jk =1 T−1∑ R−{ri} dl i jk,ri∈R,tj∈T,k∈K,l∈R−{ri},(6.3) This way we just average the distances as can be seen Fig. 6.4). Figure 6.4 Example of estimation of di jk for robot ri. In the first step we compute the predicted positions of tasks in k=2 (in green) and then, we compute all the possible distances to the predicted positions and iterate. This way we find that the distance in step k=2 of robot r1 to task t1 is the mean of the distance from 6.1 Problem Statement 101 the other tasks in their predicted positions of the previous step k=1. In k=3 d112 =d3 112 +d2 112 2,(6.4) d122 =d1 122 +d3 122 2,(6.5) d132 =d1 132 +d2 132 2,(6.6) and, the same way, in k=3we have d113 =d3 113 +d2 113 2,(6.7) d123 =d1 123 +d3 123 2,(6.8) d133 =d1 133 +d2 133 2.(6.9) In order to compute τA k , it is necessary to estimate the time that each step k∈K takes. We can compute the mean time to complete the k -th mission considering the time that robot rineeds to reach task tjas its k-th mission as ti jk =di jk vi ,ri∈R,tj∈T,k∈K,(6.10) starting from the initial position of both robots and tasks as τk=1 T·R·∑ ri∈R ∑ tj∈Tti jk +ψjk∈K.(6.11) Then, we can use τk to compute the accumulated elapsed time from the moment the allocation is performed to the beginning of task kas τA k= k−1 ∑ µ=1 τµk∈K,(6.12) and then, find the new ti jk, and use it to repeat this process until k=K. Finally, we can approximate the time it takes to accomplish a task tj in a given allocation δ,τj, by 102 Chapter 6. Predictive MRTA with dynamic tasks τ′j=∑ ri∈R ∑ k∈K t′ i jk ·δi jk,tj∈T,(6.13) with t′ i jk =τA k+ti jk +ψj,ri∈R,tj∈T,k∈K.(6.14) On the other hand, the distance traveled by the robots can be approximated as ρ′ i=∑ tj∈T ∑ k∈K di jk ·δi jk,ri∈R.(6.15) The complete algorithm to predict the positions of robots and tasks in the allocation horizon is described in Algorithm 4. Algorithm 4: di jk,ti jk,τA kand t′ i jk estimation Let k=1; Initialize the position of tasks using their current position, q′jk =qj; Compute ℓi jk using pi,q′jk,vi, and ˙qj; Compute di jk =||ℓi jk −pi||2ri∈R,tj∈T; while k≤Kdo Compute ti jk by means of (6.10); Use ti jk and ψjto compute τkusing (6.11); Compute τA k+1by means of (6.12); Estimate the position of tasks in k+1,q′jk+1, using q′jk,˙qjand the time obtained in τA k; Compute t′ i jk,ri∈R,tj∈T, by means of (6.13); Compute di j(k+1)by means of (6.3); end Taking the τ′j and ρ′ i computed using (6.13) and (6.15) the cost function can be transformed into J′=∑ ri∈R ∑ tj∈T ∑ k∈K k·(λi·di jk +ϕj·t′ i jk)·δi jk,(6.16) where k has been added to prioritize allocation δi jk with lower k so that there are no unfeasible allocations that imply a certain δi jk =1 without the corresponding variables being active, δi j(k−1)=1. 6.1 Problem Statement 103 Let us now consider the constraints in the assignment of tasks. Depending on K , we will consider two cases: 1. K·R≥T : Robots can perform all tasks during the allocation horizon. Then, we can force all the tasks to be performed once in the allocation horizon by setting the constraint: ∑ ri∈R ∑ k∈K δi jk =1,tj∈T.(6.17) Also, a robot can only perform one task for its k-th slot or stay idle by: ∑ tj∈T δi jk ≤1ri∈R,k∈K.(6.18) 2. K·R<T : Robots cannot perform all tasks in the prediction horizon. Then, robots must perform one task at each k∈K: ∑ tj∈T δi jk =1ri∈R,k∈K,(6.19) with allocated tasks being performed once at most: ∑ ri∈R ∑ k∈K δi jk ≤1,tj∈T.(6.20) Either way, it must be taken into account that a robot cannot be allocated to any task if its battery level is too low. To this end, βi jk is defined so that βi jk =1 if robot ri has enough SOC to perform task tj as its k mission and 0 otherwise. Using the prediction in (6.13), we define βi jk =(1if bi−wi·t′ i jk ≥0 −1if bi−wi·t′ i jk <0,ri∈R,tj∈T,k∈K.(6.21) and 104 Chapter 6. Predictive MRTA with dynamic tasks θik =         1if T ∑ j=1 (βi jk +1)>0 0if T ∑ j=1 (βi jk +1) = 0 ,ri∈R,k∈K,(6.22) (6.23) where θik is equal to 1 if robot ri has enough SOC to perform at least one of the tasks at k , and 0otherwise. This way, (6.19) can be relaxed as T ∑ j=1 δi jk =θik,ri∈R,k∈K,(6.24) and (6.18) can also be converted into T ∑ j=1 δi jk ≤θik,ri∈R,k∈K.(6.25) 6.1.1 Problem Formulation The problem to solve can be formulated as min δJ′,(6.26) s.t(6.17),(6.25),and (6.27) 0≤δi jk ≤1ri∈R,tj∈T,k∈K(6.28) (6.29) in the case R·K≥Tand min δJ′,(6.30) s.t(6.20),(6.24),and (6.31) 0≤δi jk ≤1ri∈R,tj∈T,k∈K(6.32) (6.33) otherwise. Note that robots will perform tasks so that the problem to be solved will be (6.26) or (6.30) depending on the number of active tasks at the time. 6.2 Predictive Multi-Robot Task Allocation Algorithm (PMRTA) 105 6.2 Predictive Multi-Robot Task Allocation Algorithm (PMRTA) In this section, the LP relaxation is described and the Predictive Receding-Horizon MRTA is detailed. Problems (6.26) and (6.30) can be transformed into min δc·δ(6.34) s.t A ·δ=b(6.35) by doing the following transformations in the cost function cT i=                           1·(λi·di11 +ϕ1·t′ i11) 2·(λi·di12 +ϕ1·t′ i12) . . . K·(λi·di1K +ϕ1·t′ i1K) 1·(λi·di21 +ϕ2·t′ i21) 2·(λi·di22 +ϕ2·t′ i22) . . . K·(λi·di2K +ϕ2·t′ i2K) . . . 1·(λi·diT1 +ϕT·t′ iT1) 2·(λi·diT2 +ϕT·t′ iT2) . . . K·(λi·diTK +ϕT·t′ iTK)                           ,(6.36) with c=c1c2··· cR,(6.37) and the following transformations in the constraints for Problem (6.26): A=               11×K01×K··· 01×K 01×K11×K··· 01×K . . .. . ..... . . 01×K01×K··· 11×K 11×K01×K··· 01×K 01×K11×K··· 01×K . . .. . ..... . . 01×K01×K··· 11×K ··· 11×K01×K··· 01×K 01×K11×K··· 01×K . . .. . ..... . . 01×K01×K··· 11×K 0T×R·K IKIK··· IK 0K0K··· 0K . . .. . ..... . . 0K0K··· 0K 0K0K··· 0K IKIK··· IK . . .. . ..... . . 0K0K··· 0K ··· 0K0K··· 0K 0K0K··· 0K . . .. . ..... . . IKIK··· IK IR·K               , b= 11×Tθ11 θ12 ··· θ1K θ21 θ22 ··· θ2K ··· θR1 θR2 ··· θRK T; (6.38) 112 Chapter 6. Predictive MRTA with dynamic tasks To generate the simulation map, we modeled a miniature version of a region in the “La Africana” thermosolar plant using built-in tools in the Gazebo simulator. The robots involved in the simulation consisted of UAVs and UGVs, which were simulated using the Parrot Bebop 2 UAV (Parrot) and Jackal UGV (Robotics) models, respectively (see Fig. 6.11). Based on the datasheet and some adjustments for the simulation, the maximum velocity of Parrot Bebop 2 UAV and Jackal UGV were set fixed as 4 m/s and 2 m/s respectively. All the models in Gazebo were then connected to the ROS environment where the algorithm was executed. (a) Close view. (b) Inclined view. (c) Upper view. Figure 6.10 Solar plant Gazebo. 6.3 Simulations and Results 113 (a) Parrot Isometric view. (b) Parrot Upper view. (c) Jackal Isometric view. (d) Jackal Upper view. Figure 6.11 Bebop Gazebo model. Then, the PMRTA algorithm was implemented and evaluated using ROS and Gazebo Simulator. GA was also used to evaluate the performance of PMRTA algorithm. In the ROS implementation, there are 5 main nodes behind the simulation. Mission Center node acts as a hub to connect the dynamics node (robot and task) and the algorithm processor node (PMRTA and GA). This node will retrieve information updates from robots and tasks, such as position, velocity, and completion status, in a hybrid fashion. Then, these information updates will be delivered to the algorithm node to recompute the allocation. Afterwards, the Mission Center node will convert the allocation into path commands for the robots. For the purpose of the Montecarlo simulation, the Problem Generator node will then randomize the initial parameter of the problems among the ranges shown in Table 6.3. This node will generate a new set of parameters when the mission is deemed complete by the Mission Center node. The complete diagram of the process is represented in Fig 6.12. Illustrative example To verify the functionality of the algorithm through simulation, a case study was conducted using 4 robots comprising 2 UAVs and 2 UGVs, along with 16 tasks moving to the east direction. The local controller for the UAVs was designed to ensure a constant altitude of 3.0 m, while the UGV local controller was designed to avoid obstacles such as the solar collector within the field. The remaining parameters for both the robots and the tasks were randomly generated and are presented in Table 6.1. For comparison, JC has been calculated for each simulation using the actual distance traveled by each robot and the actual time required for the tasks to be completed. In Table 6.2, the impacts of the allocation horizon K on the outcomes and the increase in computational load with K are evaluated. From this simulation, for K=3 and K=4 provide the best results, the same as in Section 6.3. Fig 6.13 also verifies that all tasks are 114 Chapter 6. Predictive MRTA with dynamic tasks properly performed. Figure 6.12 ROS and Gazebo architecture. Table 6.1 Illustrative example parameters. Robots pi(0) [m]vi[m/s2]λiλi(0)wi r1∈RA[5.0,2.0]4.0 2 96 0.14 r2∈RA[−2.0,−4.0]4.0 2 92 0.12 r3∈RG[−2.0,3.0]2.0 1 84 0.05 r4∈RG[3.0,−1.0]2.0 1 62 0.06 Tasks qj(0) [m]˙qj[m/s2]ϕjψj t1[−2.0,4.0] [0.11,−0.01]1 5 t2[−3.0,8.0] [0.12,−0.02]5 3 t3[−1.0,3.0] [0.09,−0.03]6 4 t4[−5.0,5.0] [0.1,−0.03]8 3 t5[3.0,6.0] [0.08,−0.01]3 2 t6[8.0,2.0] [0.09,−0.02]1 4 t7[4.0,1.0] [0.12,−0.03]2 3 t8[2.0,3.0] [0.1,−0.04]6 4 t9[1.0,−3.0] [0.08,−0.03]8 1 t10 [5.0,−2.0] [0.09,−0.06]3 2 t11 [3.0,−6.0] [0.1,−0.01]5 3 t12 [9.0,−4.0] [0.11,−0.04]9 4 t13 [−1.0,−3.0] [0.07,−0.02]3 3 t14 [−5.0,−2.0] [0.09,−0.05]6 3 t15 [−3.0,−7.0] [0.1,−0.02]4 1 t16 [−5.0,−4.0] [0.12,−0.04]1 4 6.3 Simulations and Results 115 Table 6.2 Results of the case study for variation of K. K 1 2 3 4 ... 16 JC19111.68 14122.88 8454.56 11884 ... 13980.48 riρi[m] r137.55 27.26 25.37 28.47 ... 27.91 r235.4 29.68 29.71 31.77 ... 32.14 r318.68 34.65 23.6 25.66 ... 26.5 r419.5 25.79 18.58 24.98 ... 20.88 tjτj[s] t116.54 54.22 50.95 67.7... 78.81 t251.44 40.72 22.96 36.37 ... 41.17 t358.27 38.01 7.16 26.28 ... 30.82 t443.38 27.81 20.92 20.55 ... 29.59 t528.89 13 37.21 43.08 ... 54.19 t618.06 60.41 54.06 57.84 ... 63.62 t77.39 54.77 44.94 49.16 ... 53.5 t838.25 46.3 11.35 30.53 ... 41.69 t923.03 12.31 8.94 8.45 ... 10.75 t10 26.86 6.78 36.66 40.83 ... 57.25 t11 30.42 19.39 21.83 21.95 ... 29.79 t12 67.84 47.3 29.8 11.96 ... 18.3 t13 12.26 17.02 35.48 37.59 ... 46.21 t14 41.95 14.73 10.35 9.87 ... 12.73 t15 15.43 5.54 33.89 29.7... 36.68 t16 10.64 86.74 48.32 50.35 ... 58.98 Average Computational Time 0.033 0.041 0.057 0.072 ... 0.574 Montecarlo Simulation Results For the evaluation of the algorithm, we created a set of 100 random problems. Then, the problems were solved using the proposed method with also varying K . An example of the Montecarlo simulation can be seen here. Alternatively, for greater clarity the same simulation has been represented using RVIZ, from ROS, and the results can be seen in Fig 6.13 or here. Note that the simulation can also be observed from the robots viewpoints. Table 6.3 Range of parameters. vim/s bi%wi%/s λi ri∈RA,RA=[1,4] [10,20] [0,50]%(0,1) [0,5] ri∈RG,RG=[1,4] [4,8] [0,100]%(0,0.5) [0,5] T ˙qjm/s ϕj [1,20] (0,4) (0,10) 116 Chapter 6. Predictive MRTA with dynamic tasks Figure 6.13 Snapshot of the evolution of the task operation time τj in the illustrative example Table 6.2, for K=4 . The remaining task operation time τj of each task is represented as the color bar above. Blue means that the task operation time has reached zero. This simulation can be seen here.. Note that the simulation using ROS and Gazebo also demonstrates that the optimal values of JC are obtained with K=3 and K=4 , thus confirming the results obtained from the simulation in Section 6.3. Additionally, as K increases, the cost function tends to deteriorate compared to the values for K=3 and K=4 . Moreover, in this case, the computation time generally increases with the rise in Kvalue. Figure 6.14 Mean JC , mean average computation burden and mean total computation burden of set of 100 random problems with different values of K. 6.4 Conclusions 117 6.4 Conclusions The proposed algorithm is especially suitable for MRTA problems with moving tasks, as it takes advantage of the information about the evolution of the tasks and the robot’ speeds and computes the allocation according to a certain predicted evolution. Then, it applies the first allocation to each robot and each time a task is completed it updates the information and recomputes the prediction in a receding horizon manner. The algorithm was first tested with Matlab simulations and then validated with ROS/Gazebo simulations. The results show that the proposed algorithm outperforms the greedy algorithm (proposed algorithm with K=1) since it takes into account a prediction of the future and not just the current moment and it also outperforms another well-known meta-heuristic algorithm such as the genetic algorithm. Nevertheless, the efficiency of our proposal seem to decrease if tasks move in different directions. Part III Direct Normal Irradiance Estimation 119 7 Bayesian Based Inference Estimation This chapter is based on the work of Martin et al. (2021b), where we study how the RSN can be managed to collect information for the control system in the most efficient manner. In particular, we integrate several layers of information comprising data measured by ground and air robots and fixed sensors such as pyrheliometers. Then, we use a Bayesian approach to update the probability that a certain area of the solar field is shadowed. In this way, robots are moved to positions where the information gathered is maximized. This approach is aligned with other works in the literature that use Bayesian Inference (BI) and information theory (Cover and Thomas, 2012) to manage RSN such as those of Julian et al. (2012b),Julian et al. (2013) and Cui et al. (2015), where mutual information is used in lockstep with consensus-based strategies for the control of a distributed fleet of robots acting as a RSN. The rest of this chapter is organized as it follows: In Section 7.1, the problem statement and assumptions regarding the thermosolar plant, the vehicle fleet, and the clouds are presented; in Section 7.2, the proposed algorithm for the spatial estimation is detailed; in Section 7.3, the case study where the algorithm is tested is described; in Section 7.4, the results of the simulations are shown and discussed. Conclusions and future investigation lines are given in Section 7.5. 7.1 Problem Statement Consider the PTC-CSP depicted in Chapter 2 where a grid of DNI measurement spots or cells located in the space between collectors is depicted. We assume that spatial information from previous measurements and wind vectors is available, e.g., obtained from a numeric model like Harmonie Arome in Spain (AEMET; Kalnay, 2003). Moreover, other sources of wind information such as fixed anemometers and sensors in the UAVs can be considered. Then, DNI data is obtained an converted into CF as indicated in Chapter 2. Discretized real values, CFR i j , are approached via a discrete probability distribution for 121