Indoor Visual Exploration with Multi-Rotor Aerial Robotic Vehicles

In this work, we develop a reactive algorithm for autonomous exploration of indoor, unknown environments for multiple autonomous multi-rotor robots. The novelty of our approach rests on a two-level control architecture comprised of an Artificial-Harmonic Potential Field (AHPF) for navigation and a low-level tracking controller. Owing to the AHPF properties, the field is provably safe while guaranteeing workspace exploration. At the same time, the low-level controller ensures safe tracking of the field through velocity commands to the drone’s attitude controller, which handles the challenging non-linear dynamics. This architecture leads to a robust framework for autonomous exploration, which is extended to a multi-agent approach for collaborative navigation. The integration of approximate techniques for AHPF acquisition further improves the computational complexity of the proposed solution. The control scheme and the technical results are validated through high-fidelity simulations, where all aspects, from sensing and dynamics to control, are incorporated, demonstrating the capacity of our method in successfully tackling the multi-agent exploration task.


Introduction
In recent years, the field of Robotics has gone through significant advancements in both software as well as hardware. One of the most notable subsets of the field is Unmanned Aerial Vehicles (UAVs). The rise in computing power, advances in energy storage technologies and the development of sophisticated control algorithms have resulted in UAVs being widely employed in research and industry alike. The rapid adoption of such platforms is well-motivated, owing to the versatility and recent accessibility of the platforms, which renders them suitable for a wide range of applications. In this work, we aim to utilize the capabilities of a fleet of multi-rotor autonomous aerial vehicles (drones) for the exploration of indoor spaces. While a single platform can accomplish such a task, it is evident how employing multiple agents results in a decrease in the overall exploration time, as well as overcomes possible limitations due to, e.g., limited battery life. In order to successfully accomplish this task, already existing, robust and proven technologies for visual sensing of the drone's environment are implemented. The path planning aspect will be tackled through a provably-correct custom framework based on Harmonic Artificial Potential Field (AHPFs). Finally, the control framework is completed through a low-level, custom reactive field tracker, which ensures that the robots track the field safely without sacrificing the full exploration guarantees provided by the AHPF-produced velocity field. Our framework is built around the Robot Operating System (ROS) [1], which enables sensing, communication and actuation of robotic systems. Owing to the aforementioned features, our algorithm exhibits superior performance with respect to smoothness and

State-of-the-Art on Multi-Robot Exploration
Fixed-base station multi-agent exploration has been extensively investigated in the literature through a variety of approaches [17][18][19][20][21]; however, communication constraints in such centralized frameworks have not been explicitly considered. Nevertheless, since communication constraints are a persistent aspect of real-world applications, which can hinder the efficacy and implementability of any method, there has recently been a focus on treating multi-agent-related problems under such constraints. The latter have been commonly treated through the notion of connectivity, which concerns the communication ability between robots and the fixed station. Consequently, a widely employed approach consists of maintaining the connectivity between the robots and the base station, directly or in a multi-hop manner, continuously during the execution of any related task. This constraint is especially important in the context of search and rescue applications, where valuable real-time information, such as image streaming, is pivotal [22][23][24]. Other approaches treat communication of acquired information to the base station periodically, in a discontinuous re-connection fashion, which can be implemented either synchronously [25] or asynchronously [26]. In the latter case, periodic connectivity may not be enforced as a hard constraint; for example, in the role-based distributed strategy [27], robots are free to explore the unknown environment with no communication restrictions. In order, however, to communicate relevant information, specifically selected rendezvous points are selected, where asynchronous updates of the relevant robot maps of the environment to the base station take place via communication relays. At the same time, stronger forms of asynchronous connectivity, e.g., line-of-sight communication [28], or a small set of behavioral and message exchanges between robots and dropped beacons [29], have also been investigated. Finally, communication constraints can also be defined with the goal of ensuring global connectivity only at the deployment locations of the robots through recurrent connectivity strategies. This results in enforcing recurrent connectivity upon the collection of new data by each robot [30][31][32]. The motivation behind such approaches mainly rests on the fact that robots can be disconnected for arbitrarily long periods of time while establishing a connection only when new information has been collected.

Proposed Method
Within the existing paradigm, our method aims at constructing a provably correct exploration scheme that employs the capabilities of a fleet of drones to explore an indoor environment in finite time in a safe manner. In order to take advantage of the multiple drones, no explicit communication network is assumed, and each drone can communicate with others wirelessly if and when visual identification between drones is established. This decentralized approach, while not explicitly ensuring optimality during the exploration task, results in a computationally attractive and easily implementable approach in practice, as it only necessitates an adequately equipped drone fleet.
Additionally, in order to formulate a computationally attractive method for real-time octomap building, we decouple the 3D mapping aspect from the navigation task; that is, the navigation framework relies solely on a computationally efficient 2D sector-based scheme. Conversely, the on-board visual information is employed to build a 3D representation of the workspace online. This results in a more robust navigation framework, owing to its simple workspace representation that is more attractive in the context of the on-board computing power of modern multi-rotor platforms. By guaranteeing full workspace exploration, the mapping aspect is also ensured to provide reliable results.
To summarize, our contributions are: • Development of an AHPF-based exploration algorithm for a multi-rotor platform, • Extending the above scheme to the multi-robot exploration problem, • Integrating the aforementioned exploration framework with a scheme for singleagent visual map-building of unknown workspaces, combined with an inter-agent information exchange aspect.
The goal is, therefore, to provide a computationally tractable solution for indoor exploration of multi-agent, multi-rotor platforms.

Problem Formulation
Consider a three-dimensional, simply-connected, static workspace W ∈ R 3 . In this case, the topological property of "simple connectivity" dictates that there are no floating obstacles inside the workspace. Further consider a set of drones, I = {1, 2, · · · , I}, operating within the workspace W, with p i ∈ W − ∂W, i ∈ I denoting the position of the i-th robot. The problem is thus initially formulated as follows: Assuming that the drones start from arbitrary positions within W, the goal is to explore the whole workspace in finite time, with no collisions with the workspace boundary ∂W or between the drones themselves.
In order to formulate a solution, we begin by assuming that there exists a superscribed sphere S 2 i ⊂ R , i ∈ I with radius R i for each drone, such that each drone lies entirely within the i-th sphere. Then, let W in f l denote the inflated workspace, i.e.: where R max = max{R i }, i ∈ I denotes the maximum radius of the set of superscribed spheres, and n(z) : ∂W → R 3 denotes the inwards-pointing unitary vector at the point z on the boundary of the workspace. Notice that by inflating the workspace, the problem is essentially transformed into the navigation of a point robot (this step can be easily incorporated into the proposed framework through computationally inexpensive methods). A necessary assumption that follows is that the inflation of the workspace boundary does not hinder the connectivity of the workspace; when taking the maximal radius for each drone's representation of the workspace (for compatibility purposes during map matching, as it will become apparent in the sequel, a large enough drone might not "fit" through narrow corridors of the inflated workspace. Therefore, we henceforth assume that W in f l remains fully connected. Nevertheless, this assumption is not that strong, as the existing platforms relevant to indoor applications are sufficiently small. This assumption can also be lifted if all drones are exactly similar, which is highly likely in practice. Further assume that each drone is equipped with a sensor that provides boundary information of the surrounding workspace, i.e., where B(p, r) is the disk with radius r centered at p and L(p, q) is the linear segment connecting p and q inclusively. Additionally, we denote by P i (t s , t e ), i ∈ I the i-th drone's path from time t s until t e . The explored region along the path P i , i ∈ I is thus defined as: whose boundary is denoted by ∂E( notes the free part and ∂E O (P i ) ⊆ ∂W denotes the obstacle-occupied part of the discovered boundary for the i-th drone. We additionally define an appropriate parametrization of the explored boundary, i.e., σ i :

Proposed Sub-Problems
In order to tackle the above problem, we break it down into three sub-problems, namely: 1.
Tracking Problem.
These sub-problems can now be addressed separately. Briefly, the Localization and Mapping problem essentially consists of obtaining an accurate representation of the dis-covered boundary ∂E i , i ∈ I and of the position p i , i ∈ I of each robot. Path Planning is addressed in the context of AHPFs, where a vector field suitable for navigation is constructed. If a robot follows the field exactly, then safe exploration can be guaranteed. However, a robot's dynamics, along with possible kinematic non-holonomic constraints, prohibit exact tracking of the field. Therefore, we treat the tracking problem separately through the formulation of a provably safe, reactive tracking controller.

Materials and Methods
We begin with an overview of the method. Based upon the existing attitude controllers of multirotor platforms (e.g., the open source ArduPilot software https://ardupilot.org/, accessed on 1 June 2022), we formulated a two-level controller. The latter is based on a high-level navigation field for exploration [33] and a low-level tracking controller based on the unicycle kinematic model [34]. Finally, inter-agent collision avoidance is guaranteed through the manual addition of the visible subset of agents to each agent's occupancy grid. Our method provides guarantees of safety and full-workspace exploration in finite time owing to previous works [33,35].

Multirotor Kinematics and Dynamics
In this subsection, we introduce the widely employed kinematic and dynamic models of multirotor platforms. Consider the exemplary multirotor robot, as depicted in Figure 1. Let B = e B x e B y e B z denote the body-fixed frame, whose origin coincides with the vehicle's center of mass. Furthermore, we define a fixed inertial frame I = e I x e I y e I z , as depicted in Figure 1. The translational and rotational dynamics of the vehicle are described by the Newton-Euler equations [36,37]: where p = p x p y p z T , I v = v x v y v z T denote the position and the linear velocity of the drone with respect to the inertial frame I, B v = u v w T is the linear velocity vector with respect to the body frame B, m denotes the mass, I R B denotes the rotation matrix used to perform rotations from B to I, J denotes the inertia matrix, and ω denotes the angular velocity of the robot with respect to the body frame B.
The external forces and torques that are applied to the vehicle are split into the following terms: where: denotes the drag forces and C d the drag coefficient matrix; • F g = m B R I 0 0 −g T denotes the gravitational force, where g is the gravitational constant; • F M = 0 0 T T denotes the total thrust generated by the motors; • M M = τ x τ y τ z T denotes the torque produced by the motors; • M d = C m ω ω are the drag moments with C m denoting the drag moment coefficient matrix; The total input thrust and moment applied to the drone depends heavily on the vehicle's geometrical characteristics and specifically on the number N of motors and the airframe configuration. In line with momentum theory, the thrust force T i and the drag moment τ i that are produced by the drone's propellers are assumed to be proportional to the square of the rotor's angular velocity, i.e., where i = 1, .., N and C T , C τ denote the thrust and drag coefficients correspondingly.
For the quadrotor, used in simulation scenarios, the total thrust and moments are computed by the following control allocation matrix: with l x , l y denoting the distance of each motor with respect to the center of mass.

Autopilot and On-Board Sensors
Modern multirotor platforms operate through the application of a low-level controller, which is a cascaded PID control structure consisting of an outer position loop and an inner attitude one. More specifically, the outer position loop is responsible for converting the reference position p d , velocity I v d (or B v d ) and heading ψ d of the vehicle to target orientation (roll φ d , pitch θ d , yaw ψ d ) and throttle. The inner attitude controller then translates the aforementioned orientation and throttle commands to motor Pulse Width Modulation (PWM) values. The state feedback is achieved by fusing sensor measurements, such as data by GPS, compass and IMU, using an Extended Kalman Filter. An overview of the control architecture is shown in Figure 2.

Localization and Mapping
In order to perform Localization and Mapping, robust existing SLAM tools will be employed. More specifically, we will employ the Hector SLAM method, as presented in [38]. This framework can be readily integrated through the available, open-source ROS package and provides odometry data for an estimate of a robot's position, along with a built-in Occupancy Grid Map (OCG), both of which are utilized for path planning and feedback control of the drones. Briefly, SLAM methods utilize sensor information (in Hector SLAM's case, a LIDAR sensor), which is processed in order to extract an estimate of the position of the sensor through matching sequentially acquired measurements. OCGs, on the other hand, consist of a 2D representation of a workspace and take the form of a matrix. Given an inertial frame of reference (e.g., the drone's take-off position), the surrounding (unknown) space is decomposed into squares (cells). The position (index pair) of each component of the aforementioned matrix corresponds to a cell inside the workspace, while the value of each component represents the probability that the respective cell is occupied by an obstacle, given the currently available information to the drone. Therefore, having such a representation, obtaining the sensed boundary ∂E can be directly tackled through well-known boundary extraction imaging techniques [39,40], applied over the OCG.

Velocity Field
In this section, we treat the path planning problem for a point-robot (as discussed in Section 3). In order to achieve provably complete path planning for each robot, we employ the AHPF-based method presented in [33]. Briefly, the method designs a reactive vector field for the single integrator dynamics, i.e., given the problem formulation of Section 3, along with the following dynamics: the authors in [33] provide a reactive vector field, which is modified herein to fit the multi-agent framework: where K u ∈ R is a constant gain, the set: for the drone i denotes all other drones that lie inside its sensor radius, S R 1 : R → R is a scaling function to ensure safety and will be defined in the sequel, and finally, d(p i ; E i ; p i I ) is a distance function: which outputs the minimum distance between the already sensed boundary by the i-th drone and its distance from the closest of the other drones. The quantity D i min > 0 is a parameter that determines the least permissible distance between two drones and is chosen as D i min < R i , i ∈ I. Furthermore, we define the function: The velocity control law (13) is further derived from the gradient of the reference potential φ(p;k) : W → R, which is determined by a set of parametersk, which will be defined in the following sub-section. The potential φ i , i ∈ I of the i-th drone is designed to satisfy the Laplace equation: with Neumann boundary conditions: where σ −1 (q) is the inverse function of the parametrization of the boundary introduced in Section 3. Effectively, the function k σ −1 (q), t dictates whether a point on the boundary q ∈ ∂E is attractive (k < 0) or repulsive (k > 0). In order for the Neumann boundary problem to admit a solution, the boundary conditions should satisfy the compatibility constraint: We, henceforth, refer to such a set of boundary conditions as compatible. The reason why we employ such a potential for navigation is that such solutions lack local equilibria inside the domain over which (17) is solved. This property holds according to the minimummaximum principle, which dictates that solutions to the Laplace equation exhibit minimal and maximal values only at the boundary of their domain. The boundary condition (18) further ensures that the resulting gradient vector field is safe, i.e., the resulting velocities point inwards at the boundary.

Boundary Discretization
It is evident that the Laplace problem (17), along with the boundary conditions (18), necessitates the definition of the function k(·), which would result in a problem with an infinite-dimensional boundary condition, as (18) is defined over a continuous set. In order to handle the boundary conditions numerically, we employ the following parametrization of k over a finite number of control pointsk, i.e., where each k j i , j ∈ {1, · · · , N C p } corresponds to points q j placed on the observed boundary. For any boundary point q j ,k j i corresponds to the exact value of the boundary condition (18) for the i-th drone, while the same value for any point q = q j , k σ −1 i (q), t,k i is chosen as a linear combination of nearby control point values. The discretization scheme is illustrated in Figure 3. In addition, we equip the specified boundary conditionsk with appropriately designed adaptive laws, such that no additional (locally) stable equilibria appear in int(∂E).
These values are tuned according to the adaptive law:k (the index i corresponding to each drone is omitted for the sake of clarity) where c, µ, b e ,k t are provided in [33] and are omitted here for the sake of brevity. Briefly, the valuesk t operate as reference values for the boundary conditions. It is evident that over time, (21) results ink converging to the reference values. Additionally, c renders the first term of (21) zero when the robot is near a critical point of the field φ while b e ensures that the trajectories of the robot remain safe for all time. Finally, µ is employed to modify the convergence rate.

Fast Multipole Boundary Elements
In order to obtain a solution to the Laplace Equation (17) given the boundary conditions (18) dictated by the function k of the preceding sub-section, we employ a Boundary Element Method (BEM) and, more specifically, the Fast Multipole Boundary Elements Method (FBEM), as employed in [33]. Briefly, in FMBEM, the computational cost of the conventional BEM (which is O(n 2 ), where n is the number of boundary elements used for approximating the domain's boundary) is significantly reduced through the computation of approximate solutions for the corresponding boundary value problem, with a specified error. In order to accomplish this, a quadtree decomposition is employed in order to obtain a hierarchical subdivision of the elements that describe the boundary of the solution domain. Subsequently, multipole expansion is employed to approximate the interaction among elements (with respect to the Laplace equation). Therefore, an analytical computation of every pairwise element interaction is avoided, at the expense, however, of an approximation error, which is nevertheless bounded. This process results in a reduction in the computational complexity down to O(n). Finally, in FMBEM, there is no need to store the dense and non-symmetric matrix that is part of the conventional BEM if an iterative solver is employed, e.g., GMRES [41], which consequently also results in O(n) memory requirements.

A Brief Discussion about the Algorithm
In this subsection, we briefly explain how the FBEM solution to the problems (17) and (18) is employed along with (21) to provides a safe navigation field suitable for exploration. As in [33], a solution at time T > 0 to the problems (17) and (18) is obtained by applying the FBEM method, along with the boundary conditions that stem from the parameters' values at that time, i.e.,k(T). Once such a solution is obtained, the linear dependence of the potential on these parameters, along with the adaptive control law (21), provide a safe field for exploration, as proven by the technical results in [33]. As each drone explores the workspace, the Laplace problem is updated according to newly-discovered parts of the boundary, and successive solutions are obtained. In the following sub-section, we provide the technical results of [33] that prove the asserted claims.

Technical Results
In this subsection, we include some of the technical results presented in [33] for completeness while omitting the relevant proofs for brevity. For details, we direct the reader to [33].

Proposition 1.
Assuming thatk is compatible at t = 0, the adaptive law (21) guarantees thatk will remain compatible for all time. φ(p,k) is unsafe, the adaptive law (21) guarantees that it will become safe in a finite amount of time.

Proposition 3.
All equilibria of the dynamic system: located in int(E) are unstable. (8) in [33]). The proposed control protocol (13) and (21) ensures that the trajectories ofṗ i = u i , i ∈ I are safe.

Proposition 4 (Proposition
The above propositions ensure that the proposed velocity field has the necessary properties that will result in safe exploration of the whole workspace for each drone in finite time. Additionally, in (13), we propose a modification with respect to [33], which includes the distances between drones that lie inside the sensor radius of each other and results in safety between drones. In order to guarantee this additional type of safety, an additional feature is proposed; namely, within the OCG of each drone, we place a "fictional" obstacle, equivalent to the superscribing sphere S 2 i that was discussed in Section 3.

Tracking Controller
While we have introduced a provably complete navigation field for each drone, we need to ensure that, despite the inability of a drone to follow the single integrator dynamics exactly (12), safety is not violated. We assume that there exists an attitude controller over the drone dynamics that takes reference velocities as inputs and results in an adequate response by the drone. Such controllers are in widespread use in the industry and are standard in any available drone platform [42,43]. Additionally, we need to ensure that the drone's visual system is also pointing towards the direction of the motion, which, besides ensuring proper sensing, results in a better response with respect to the commanded velocities, as the drone is not commanded to alter its direction of motion in a sharp, non-smooth fashion. If the drone's low-level tracking controller was commanded to change its direction in such a fashion, this would result in a delayed response (owing to the drone's dynamics and its momentum), which might render its motion unsafe. To this end, we first introduce the yaw angle of the drone θ ∈ S 1 (where S 1 is the set of 1-dimensional rotations), which can be obtained with respect to an inertial frame of reference {I} as θ = atan2( I p x , I p y ), where [ I p x , I p y ] denotes the x and y axes of the drone's own frame of reference expressed in {I}. Then, we propose the following reactive control law: where α is the relative error angle, i.e., where θ ref is the reference angle of the underlying velocity field, i.e., We prove how this control law stabilizes the reference angle at α = 0 ⇔ θ = θ ref : Proposition 5. The controller (23) results in zero relative error angle over time for a constant reference angle θ ref .
However, it is evident thatθ ref = 0 as the drone moves inside the workspace, which changes the vector of the reference velocity and, therefore, its reference angle as well. To remedy this, we propose the following velocity controller: where γ ∈ (0, 1), σ p , σ u : R + → [0, 1] are smooth and monotonic bump functions such that σ p (0) = 0, σ p (∞) = 1. The quantityê i n denotes the unitary vector along the direction of the i-th drone's heading (forward-looking direction),n(p) denotes the unitary, inwardspointing vector at the closest point of boundary to the drone's position p, and ·, · denotes the classical inner product between two vectors, while u(p) denotes the reference velocity field presented in Section 4.4. More specifically, we define the normal vector at the i-th drone point p i :n (p i ) = n(z(p i )), z(p i ) = arg min where evidently n(z), z ∈ ∂E i ∪ p i I denotes the normal vector to the boundary of the so-far explored region or the closest of the rest of the robots (if it lies within the i-th robot's sensor radius). The first term in (28) becomes zero as the drone approaches the boundary of another drone (where d(p; E i ; p i I ) → 0). The second term approaches zero as: • The vectorsê n , u(p) point to different directions, and • The vectorsê n ,n point to different directions.
The first condition means that the velocity of the drone decreases if the robot is misaligned with respect to the velocity field, while the second condition decreases the velocity if the drone is pointing outwards at the boundary. Additionally, this implies that if the functions σ p , σ u , along with the tuning parameters γ, p , u , are chosen appropriately, they provide an adequate bound toθ, which renders the controller (23) asymptotically stable forθ ref = 0. This results in safe and successful navigation, regardless of the drone's dynamics, if the drone's own attitude controller can perform adequately. The efficacy of the above-proposed control protocol in successfully completing the multi-drone exploration task is demonstrated in rigorous, high-fidelity software-in-the-loop (SITL) simulations.

Octopmap Building
Since the navigation aspect is fully addressed through a two-dimensional OCG framework, the visual information is employed independently in order to build a threedimensional representation of the indoor environment. This is achieved through integrating already existing, robust tools for 3D mapping. The independence of the navigation and 3D mapping modules results in a computationally more attractive scheme. We demonstrate the efficacy of the proposed method in building a complete map in the results section of our work in a realistic, synthetic simulation environment.

Drone Communication and Map Merging
In order to exploit the full potential of the multi-drone exploration framework, a method for establishing communication between the drones to exchange the robots' maps needs to be formulated. Subsequently, each robot should process the exchanged map information and merge the latter with its own map in order to take advantage of the exploration of the rest of the drone fleet. To accomplish this, we formulate a decentralized, aperiodic communication strategy. In order to provide a simple solution that negates the need for a centralized station, we assume that during the exploration process, a number of C ∈ N drones, denoted by {i 1 , i 2 , · · · , i C } C ⊂ I, lie in the visual field of one or more of the aforementioned drones. We formulate the communication scheme for one drone, denoted by i ∈ I, and employ the same strategy for any other drone. In order to verify the presence of the set of drones C i ⊂ I − {i} for the i-th drone within its vicinity, we use visual means. This choice further motivates the Lyapunov-based controller for each drone's heading angle θ (23), as it is necessary that each drone's visual system is appropriately oriented towards its direction of motion in order to possibly detect any other neighboring drone. We further assume that the presence of two drones within the field of view of any of the latter implies that wireless communication is feasible. Having successfully satisfied the above conditions, the i-th drone can send its formed OCG map towards every drone j ∈ C i in its vicinity.
A pivotal aspect to the success of the scheme is the merging, by the i-th drone, of a composed map from the set of neighboring drones C i . There exist several methods in the literature that achieve this merging, from [44], where SLAM was extended for a set of multiple robots based on an Extended Kalman Filter (EKF), to [45], which employed a particle filter to achieve multi-robot SLAM. In another approach [46], the author proposes a hybrid method for the formation of a maximum likelihood map along with a Monte-Carlo position estimator. Furthermore, there are off-the-shelf ROS packages available for map-merging, e.g., [47]. In this work, we employ this tool to achieve the map merging for each drone.

Results
In order to demonstrate the efficacy of the proposed scheme, we have constructed a synthetic simulation scenario where two drones operate within an a priori unknown indoor environment, comprising of planar vertical walls. The environment is depicted in Figure 4. As the drones navigate autonomously, they gather visual data and build an Octomap 3D representation of their environment. Successive snapshots during the exploration process are depicted in Figure 5. We note that both Octomaps, one formed by each drone, are merged and presented in the respective snapshots for illustration purposes. The final Octomap, along with the merged OCG maps after the exploration process has been completed, are presented in Figure 6. Furthermore, the trajectories of the two drones overlayed on top of the final OCG map are depicted in Figure 7. Finally, a video of the experiment can be accessed through the following hyperlink: https://youtu.be/yDP_Zz-P0lE (accessed on 11 June 2022).

Limitations
The main limitations of this work are related to sensing and tuning. Most multi-rotor platforms are not equipped with an omni-directional Lidar, but rather their sensing is limited over a smaller subset of the full 360 • range. This specification might result in the formal assumptions of the method being violated, owing to the drone "exiting" the sensed area over which the AHPF is defined, evidently due to the drone's blindspots. This issue can be addressed if, prior to initiating exploration, the drone performs a full rotation about the vertical axis of its body frame of reference, thus forming an appropriate "sensed region" around itself and negating the probability of exiting the aforementioned area.
The second limitation is due to the use of Lyapunov/PID tracking controllers. Since these controllers are model-agnostic, some tuning is required to ensure proper tracking of the underlying AHPF. Nevertheless, tuning the respective parameters did not prove to be difficult in practice, especially since a conservative choice of maximum linear velocity can easily result in safe navigation, however, sacrificing total exploration time. An interesting future direction rests on exploring tuning methods for the parameters of the herein proposed controllers.
Another envisioned limitation-although not observed in practice during the evaluation of the proposed scheme-lies in the computational complexity of the method in large workspaces. While the FBEM is advantageous insofar as it scales linearly with respect to the size of the boundary (in contrast to scaling proportionally to the area of the workspace), very large workspaces could still pose an issue as the whole workspace is considered during the solution to the Laplace problem. This could evidently be sidestepped, as parts of the explored workspace could be discarded. This is another interesting future research area where we intend to formulate a polygonal convex decomposition method for reducing the computational load of the proposed scheme in large workspaces.
A final limitation stems from the limited battery life of drones, which is crucial in the context of exploration, as the size of the workspace is not known a priori. This could be addressed through existing approaches, such as [48].

Conclusions
In this work, we have formulated a method for solving the cooperative exploration task through a multirotor-based scheme. The proposed framework relies on existing robust tools for mapping and localization, while its novelty lies in the application of a Harmonic Potential-based exploration algorithm. Through the technical results, along with the in-simulation validation of the proposed control scheme, our method effectively tackles the multi-agent exploration problem with multi-rotor aerial robots. We, therefore, conclude that the herein posed research problem was successfully solved. Furthermore, our method exhibits some significant advantages. Since the proposed scheme employs the FBEM method, transforming the problem into a boundary problem on a 1D curve, it scales favorably with respect to the size of the workspace (O(n) complexity). This is in contrast to other popular methods, such as other popular planners for navigation in unknown workspaces, such as A [49], which, at best, scales polynomially or, at worst, exponentially with respect to the size of the workspace. Another advantage of the method is the two-layer approach through a tracking controller over a navigation field. This results in a more robust scheme with respect to the drone's dynamic behavior, as the internal attitude controller of the multi-rotor platform is employed to tackle the stabilization of its highly non-linear dynamics. Furthermore, our method is modular, which means that it can be easily employed in several robotic platforms (e.g., mobile robots) with minor modifications. This results in a robust and computationally efficient framework for cooperative exploration of indoor areas, as demonstrated in a high-fidelity realistic simulation scenario.
With regards to future work, we intend to improve several aspects of the algorithm. First of all, the parameters of the proposed tracking controller are tuned by hand. It would, however, be possible to formulate tuning laws based on some optimality criterion, e.g., exploration time, explored area coverage, etc., in a multi-agent, game-theoretic approach.
Additionally, since computing the field for the whole workspace is unnecessary in several cases, thus adding to the computational complexity of the method with no apparent advantages, we intend on formulating an algorithm for disregarding subsets of the explored workspace that do not "contribute" to the discovery of new parts of the workspace, based on polygonal decomposition. This will effectively bound the number of parameters as in [50].