Coverage Path Planning with Budget Constraints for Multiple Unmanned Ground Vehicles

Vu Phi Tran, Asanka Perera, Matthew A. Garratt, , Kathryn Kasmarik, , Sreenatha Anavatti

Abstract

This paper proposes an innovative approach to coverage path planning and obstacle avoidance for multiple Unmanned Ground Vehicles (UGVs) in a changing environment, taking into account constraints on the time, path length, number of UGVs and obstacles. Our approach leverages deformable virtual leader-follower formations to enable UGVs to adapt their formation based on both planned and real-time sensor data. A hierarchical block algorithm is employed to identify areas in the environment where UGV formations can spread out to meet time and budget constraints. Additionally, we introduce a novel control scheme that allows each UGV to generate a local steering force to dodge any static and mobile obstacles based on the closest safe angle. Results from simulations and real UGV experiments demonstrate that our approach achieves a higher coverage percentage than rule-based and reactive swarming approaches without planning. Our approach offers a promising solution for efficient coverage path planning and obstacle avoidance in complex environments with multiple UGVs.

Index Terms:

Coverage Path Planning, Spanning Tree Coverage, Optimisation Technique, Formation Control, Obstacle Avoidance, Autonomous Vehicles.

I Introduction

Intelligent transportation is gaining popularity with an increase in the number of practical applications [1]. This is particularly true for autonomous systems, such as unmanned ground vehicles (UGVs). UGVs have different applications, including coverage of a given area. In particular, when a number of UGVs/agents are employed to cover a given area, the control systems need to be intelligent to achieve the mission while overcoming the obstacles, both static as well as dynamic. The coverage path planning problem is a task wherein a UGV or UGVs, possessing a complete geometric description of the area of interest, generates an efficient coverage path to visit every point in a given area while avoiding all possible obstacles [2]. Various technological developments and advancements in sensor technology, navigational, communication, and computational systems have facilitated the rapid growth in the use of coverage path planning (CPP) methods to assist UGVs in performing many specific applications, ranging from humanitarian missions such as surveillance, search and rescue tasks, to military operations such as surveillance [3], environmental monitoring [4], and civilian applications such as area cleaning, seeding or harvesting [5], and mapping and model reconstruction [6].

In recent years, the literature has discussed several approaches for coverage by a single vehicle [7]. However, real-world factors such as battery capacity or sensor payload restrictions [8] may limit the ability of a single agent to meet an operational time limit. Compared to single-vehicle CPP, a group of multiple vehicles may solve a coverage task more rapidly due to its larger footprint. Yet, to exploit the capacity of a multi-vehicle team, novel algorithms are required to determine the route for each vehicle when they can spread out and take on a narrow formation.

Motivated by the aforementioned observations, the contributions of the present paper are:

•

A novel problem definition for non-backtracking coverage path planning with budget constraints assuming the use of a multi-UGV team in cluttered and uncertain environments.
•

A novel algorithm is presented for solving the problem, which utilizes a hierarchical block approach to decompose a given map into appropriate cell sizes. This allows us to exploit flexible multi-UGV formations to meet multiple budget constraints.
•

A distributed virtual leader-follower formation control strategy including automatic role assignment in the formation and obstacle avoidance.
•

A comprehensive comparative study in both simulated and real-world settings to confirm the viability of our approach. Our approach outperforms existing methods in terms of maximum coverage percentage, time to achieve coverage and computational complexity.

The virtual leader-follower control approach taken in this paper is based on the virtual spring system [9]. The virtual leader-follower approach offers several advantages over the traditional leader-follower and swarm-based methods in coverage path planning. Firstly, it eliminates the need for a physical leader, which can be costly and risky to implement. Secondly, it allows for the efficient coordination of multiple followers without the risk of collisions or formation breakage, which is a common issue with swarm-based approaches. Thirdly, a virtual leader can easily adapt to changes in the environment, dynamically adjust the path plan, and provide more accurate and reliable instructions to the followers. This is in contrast to the real leader-follower system, which may not be able to respond quickly enough to changes in the environment [10, 11]. Finally, the virtual leader-follower approach offers more flexibility and scalability, enabling the coordination of a large number of followers without requiring additional resources.

The remainder of this paper is organised as follows. Section II discusses related work from the literature. Section III states our coverage path planning problem definition and describes our approach to solving this problem. Our approach has components for path planning and prediction of how long it will take to follow a path in formation. Section IV presents a series of experiments with each of these components, first in simulations then on real UGVs in an outdoor setting. We offer conclusions and directions for future work in Section V.

II Literature Review

Various techniques for area coverage exist, including Voronoi-based [12], graph-based [13], next best view [14], frontier-based [15], and spanning tree coverage methods [16]. While some techniques provide incomplete coverage, others like spanning tree coverage guarantee complete coverage by generating a non-overlapping path along the spanning tree. In this study, we employ the spanning tree method for multi-unmanned ground vehicle (UGV) settings where several UGVs are used to cover an area.

Developing Multi-Cell Path Planning (MCPP) strategies for real-world scenarios is a challenging task that requires considering additional factors and requirements. Previous studies on MCPP have mainly focused on obstacle-free spaces [17, 18] or single cover types [19], which may not apply to many real-world scenarios. Only a block data structured method [20] has been developed so far to obtain a comprehensive and safe area coverage path for a team of vehicles, using the concepts of a contour map, connected graph, and spanning tree. However, this method may not scale well to large-scale multi-vehicle systems due to the tiny cells generated around obstacle boundaries. Additionally, no existing methods focus on dynamic environments with non-stationary obstacles, which are more common in practice. Therefore, detecting and avoiding moving obstacles correctly and efficiently is crucial.

Most MCPP algorithms rely on centralized control, which can lead to communication overhead and coordination difficulties as the number of agents increases [21]. To address this, a distributed approach can be used to improve coordination and communication among vehicles, generate more efficient coverage paths, and better adapt to changes in the environment. Additionally, a distributed approach provides greater robustness to failures or disturbances (e.g., the loss of a vehicle, changes in the environment, or communication failures) by allowing the system to reconfigure dynamically [22].

Furthermore, most MCPP algorithms, such as multi-robot forest coverage [23] and spiral spanning tree coverage [24], assume a desire for 100% area coverage without recharging or refilling robots [25]. However, in real-world applications, many physical constraints need to be considered, such as limited exploration time, travelled path length, restricted access to some areas, and different types of coverage required. Therefore, the state-of-the-art MCPP strategies compute plans that lead to imperfect coverage but are considerably more energy/time-efficient [26]. Additionally, complete coverage may not even be feasible when areas to be covered have impeded or hidden components, further highlighting the need for partial coverage [27].

Several solutions have applied multi-objective evolutionary techniques in path planning to balance coverage and energy/time [26, 28]. However, such multi-objective approaches cannot provide a satisfactory solution when multiple objectives are considered, as in many real-world scenarios. To the best of our knowledge, optimizing the grid cell size, a key parameter in the MCPP problem with constraints, while meeting all specific goals and limitations, has not been extensively examined. Additionally, for large-scale workspaces, multi-objective evolutionary optimization techniques are often time-consuming due to offline training. Therefore, the presented methods are not appropriate for a dynamic setting, where real-time decisions must be made [29]. A simpler alternative with very low computational complexity must be taken into account.

Moreover, controlling a flexible leader-follower formation to track a pre-defined path is a relatively unexplored approach compared to allocating a specific area to each vehicle. However, this approach has significant benefits which we wish to exploit. By staying close together in a re-configurable formation, the robots can communicate and use visual relative positioning, which can be crucial in scenarios where GPS signals are degraded or denied. The leader-follower approach also allows the vehicles to adapt to changes in the environment or mission requirements, which is not possible with fixed area allocation. This approach also enables the formation to be more robust to failures or disturbances, as the leader can quickly adjust the formation to account for the loss of a vehicle, making the system more resilient [30].

Overall, the MCPP problem is a complex optimization problem that requires consideration of multiple factors and constraints, including area coverage, energy/time efficiency, communication overhead, robustness, obstacles, and uncertainty. While there have been many advancements in this field, there is still much work to be done to develop practical and effective MCPP strategies that can be applied to a wide range of real-world scenarios. In the next section, we formulate the problem addressed in this paper and present our solution.

III Problem Definition and Algorithm

In this section, we begin in Part A by describing the MCPP problem under the physical limits we address in this paper. We provide our algorithmic solution in Part B, followed by a complexity analysis in Part C.

III-A Problem Formulation

The simple unicycle model captures the trade-off between linear and angular velocity control for small unmanned ground vehicles (UGVs) [31] and forms the basis of the UGV simulation used in this paper. We denote $q$ as the UGV state comprising its position and heading, whilst $u$ is the velocity command vector:

q=[x~{}y~{}\theta]^{T}\in\mathbb{R}^{3},u=[v~{}\omega]^{T}\in\mathbb{R}^{2}.

(1)

We will refer to two components of the state $q$ : $q_{p}\triangleq[x~{}y]^{T}\in\mathbb{R}^{2}$ , the position, and $q_{\theta}\triangleq\theta\in\mathbb{R}$ the yaw angle. The coordinates $x$ and $y$ capture the position of the center of gravity of the vehicle while $\theta$ is the angle of the wheel with respect to the $x$ -axis of the reference coordinate frame. The control variable $v$ is the forward velocity of the vehicle (along its body fixed frame $x$ -axis) while $\omega$ is the rate of change of the yaw angle.

The following equations capture the unicycle model with the known initial conditions $x_{0}$ , $y_{0}$ and $\theta_{0}$ :

\begin{split}\dot{x}&=v~{}cos(\theta),x(0)=x_{0},\\ \dot{y}&=v~{}cos(\theta),y(0)=y_{0},\\ \dot{\theta}&=\omega,\theta(0)=\theta_{0}.\\ \end{split}

(2)

We assume the environment’s shape and obstacles are known. Consider an environment $\mathcal{C}$ to be explored by a formation of $n_{q}$ vehicles, $\mathcal{Q}=\{q_{1},...,q_{n_{q}}\}$ . We choose a tessellated environment of grid cells of width $C_{W}$ and height $C_{H}$ such that $\mathcal{C}=\{c_{ij}|i=1,\dots,\mathcal{C}_{W},j=1,\dots,\mathcal{C}_{H}\}$ . Given a discretised map $\mathcal{C}$ with a set of known obstacles $\mathcal{O}$ , a leader-follower formation with initial configuration $Q^{0}\in\mathcal{C}_{free}$ and $\mathcal{C}_{free}\triangleq\mathcal{C}\setminus\mathcal{O}$ , a time budget of $T_{max}$ steps, (and/or a UGV path length budget of $L_{max}$ ) and a set of observed cells at time $k$ , $\sigma^{k}:=\sigma^{k-1}\cup\zeta^{k},$ where $\zeta^{k}=\cup_{i=1}^{n_{q}}(\zeta_{i}^{k})$ is the set of new cells covered by all UGVs’ total observation area at the time $k$ . Calculate a plan $\mathcal{P}^{T}:=Q^{0},...,Q^{k},...,Q^{T},\forall T\leq T_{max}$ that does not cross any obstacles $(Q^{k}\in\mathcal{C}_{free},\forall k)$ , and maximizes the coverage percentage $CP$ with respect to cell size $CS$ . In other words, the MCCP problem with constraints requires us to maximise:

CP=\frac{\sigma^{T}}{C_{free}+OB}\times 100\%

(3)

subject to:

\begin{split}&\sum_{k=0}^{T}\lVert Q_{p}^{k}-Q_{p}^{k+1}\rVert_{2}\leq L_{max},\\ &\lVert Q_{p}^{k}-Q_{p}^{k+1}\rVert_{2}\leq v_{max},\forall~{}k,and\\ &T\leq T_{max}.\\ \end{split}

(4)

where $Q_{p}^{k}\in\mathbb{R}^{2}$ represents the coordinates of the centre of the formation. Additionally, $OB$ is the areas occupied by the obstacle cells, $L$ and $T$ denote the actual path length and coverage time, respectively. $v_{max}$ indicates the maximum velocity of the mobile vehicle. The nomenclature for the problem and our algorithm is summarised in Table I.

TABLE I: Nomenclature

Symbols	Parameters
$L_{max}$	Path length budget
$T_{max}$	Coverage time budget
$CP_{max}$	Maximum coverage percentage
L	Actual path length
T	Actual coverage time
CP	Actual coverage percentage
MST	Minimum spanning tree
$\hat{L}$	Predicted path length
$\hat{T}$	Predicted coverage time
$\hat{CP}$	Predicted coverage percentage
$n_{q}$	Number of mobile vehicles
$\theta$	UGV heading
$q_{p}$	UGV position in 2D
$Q_{p}$	Centre position of the formation
CS	Grid cell size
BS	Block size
OB	Obstacle cells located in the UGVs’ total observation range
$\zeta^{k}$	New cells covered by all UGVs’ total observation area
$\sigma^{k}$	Total observed cells at time $k$
$l_{o}$	Desired natural length
$l_{a}$	Actual length
$F_{ij}$	Virtual spring force vector between UGV $i$ and UGV $j$
$D_{c_{ij},c_{i^{{}^{\prime}}j^{{}^{\prime}}}}$	Connection direction from cell $c_{ij}$ to $c_{i^{{}^{\prime}}j^{{}^{\prime}}}$
$F_{g}$	Goal following force vector
$\omega_{g}$	Target force weight
$\Delta$	Blocked angles

III-B Algorithm for Coverage Path Planning using Flexible Formation Control

This section describes our algorithm for MCPP in a known environment. We begin with an informal description supported by diagrams here, then provide algorithmic details in the following sub-sections. The algorithm is designed to produce a coverage path that maximises the area of the environment that will be visited while meeting a given time budget. We use a binary and linear search to find the smallest grid cell size (CS) that will produce a coverage path in the presence of obstacles as per (3). The algorithm proceeds as follows:

1.

The algorithm starts by choosing a grid cell size and overlaying it on the given map. For example, grid lines can be seen in grey in Fig. 1.
2.

Next, a scanline algorithm is used to convert the environment geometry onto the grid as filled or free cells. The blue cells in Fig. 1 show (filled) obstacles, while yellow cells represent (free) open space.
3.

The open space is then decomposed into variable-sized square blocks that follow the grid lines. These blocks are colored red, aqua, and purple in Fig. 1.

(a)

(b)

Figure 1: (a) demonstrates the coverage path obtained with a fine grid of 1.4m, resulting in a longer path; (b) presents an example with a coarse grid of 2.9m, leading to a shorter coverage path. The UGVs navigate through waypoints positioned at the bends of the path. The extent of the environment covered by the UGVs depends on the grid cell size due to the grid discretization. Blue polygons depict the obstacles, while purple cells indicate areas occupied by obstacles. Yellow cells represent the traversal areas, and purple cells with yellow lines depict the boundaries of obstacle areas. The small circles indicate sharp turns in the path. The connected green lines illustrate the robot path planned using the minimum spanning tree (MST) technique.
4.

The centre points of these blocks are joined by a minimum spanning tree, shown in black in Fig. 1.
5.

Afterward, a coverage path is constructed that circumnavigates the spanning tree. This is shown as a dashed green line in Fig. 1.
6.

Using the coverage path, the algorithm predicts how long coverage will take. Based on this predicted time and maximum path length, the binary search either outputs the final grid cell size or continues to dot point one above until a grid cell size is determined that will meet the budget. If the cell size increases above the observation range of the UGV, the algorithm will exit without a solution.
7.

Once the cell size is chosen, waypoints are generated at each bend on the coverage path. These are shown as black circles in Fig. 1 and contain data about the size of the block they sit in.
8.

The waypoints are then sent, in order, to a group of UGVs, which use the block size data to assume a formation that will enable them to pass through the block in a way that covers all the cells within it. The group of UGVs adapt their formation in response to both planned and real-time data from their sensors. Example formations are shown in Fig. 5.

Fig. 1(a) has a fine grid, while Fig. 1(b) has a coarse grid. The effect of the grid resolution is twofold. First, the obstacle boundaries (shown in dark purple cells) become coarser, blocking off (or de-prioritising) parts of the environment. For example, in Fig. 1 the area behind the large obstacle is blocked off. This, in turn, has the effect that the spanning tree becomes smaller, and the coverage path is shorter. These two effects contribute to the UGVs being able to meet a lower time budget by trading-off coverage percent. Further, Figure 2 shows the flowchart of the algorithm, which provides a visual representation of the steps described in the algorithm list. We now consider the algorithm in detail.

Refer to caption — Figure 2: Flow chart of the whole strategy. See page 3 for details.

III-B1 Interpreting Obstacle Geometry as a Grid

We assume obstacles are input to the system as a set of vertices. A scanline algorithm [32] is used to classify grid cells as being either an obstacle or free space. The scanline algorithm scans through the grid horizontally. Locations of the intersection between the scan lines and the obstacle edges are detected. Each obstacle is then discretized into grid cells at all intersection points and between pairs of intersection points.

III-B2 Creating Variable Sized Blocks

Once the accessible areas are realized, a block-building algorithm is applied to group adjacent free grid cells into variable-sized blocks. The largest block size is chosen to accommodate the size of the group of UGVs when they are spread out. The remaining block sizes are calculated by progressively halving the largest size.

We assume several window-based scans are performed from left to right and from bottom to top using each block size in turn. The largest block size is selected first, then gradually reduced to 1. If all grid cells in the scan window are free and do not belong to another block, they will be merged into a block and labelled with the block identifier. Otherwise, the scan window will skip and continue to the next position. After scanning all block sizes, a list of blocks covering the entire map is obtained. The block-building algorithm is summarised in Algorithm 1 of the Supplementary Material section.

Next, the connectivity between the blocks needs to be established. Every grid cell always has four connected neighbours (top, bottom, left, right). See Fig. 3 for an example. The white and black cells are free and obstacle cells, respectively. The black lines are the edges of the spanning tree. The green dash lines are the generated paths. Blue dash lines and yellow dots are the boundary lines and the center of block parts. If the two adjacent cells belong to two different blocks, these two blocks will be marked as connected. Meanwhile, the connection direction of the two blocks is equal to the connection direction of the two cells. However, there is also no connecting edge between the center block and the bottom left block in the MST (black line) because the top left block is closer to the bottom left block than the center block. The connection direction $d\in\{TOP,LEFT,BOTTOM,RIGHT\}$ , from cell $c_{ij}$ to cell $c_{i^{{}^{\prime}}j^{{}^{\prime}}}$ $d_{c_{ij},c_{i^{{}^{\prime}}j^{{}^{\prime}}}}$ is derived as:

d_{c_{ij},c_{i^{{}^{\prime}}j^{{}^{\prime}}}}=\begin{cases}LEFT~{},~{}[i-i^{{}^{\prime}},j-j^{{}^{\prime}}]=[-1,0],\\ RIGHT~{},~{}[i-i^{{}^{\prime}},j-j^{{}^{\prime}}]=[1,0],\\ BOTTOM~{},~{}[i-i^{{}^{\prime}},j-j^{{}^{\prime}}]=[0,-1],\\ TOP~{},~{}[i-i^{{}^{\prime}},j-j^{{}^{\prime}}]=[0,1].\\ \end{cases}

(5)

The block to which grid cell $c$ belongs is $b_{c}$ . Assume blocks of adjacent cells $c_{ij}$ and $c_{i^{{}^{\prime}}j^{{}^{\prime}}}$ are different. The direction from the block of cell $c_{ij}$ to the block of cell $c_{i^{{}^{\prime}}j^{{}^{\prime}}}$ is given by:

d_{b_{c_{ij}},b_{c_{i^{{}^{\prime}}j^{{}^{\prime}}}}}=d_{c_{ij},c_{i^{{}^{\prime}}j^{{}^{\prime}}}}

(6)

After building a graph with the connection between blocks, the minimum spanning tree algorithm will be used to find a spanning tree. Connections that do not belong to the spanning tree will be removed.

III-B3 Minimum Spanning Tree

A spanning tree is a subset of an undirected graph $G(V,E)$ that connects all the vertices $V$ of the graph with a minimum number of edges $E$ . A spanning tree cannot contain cycles or disconnected nodes. However, a connected and undirected graph $G$ may be connected with more than one spanning tree since every vertex can be connected from many directions. The cost of the spanning tree is the sum of edge weights in the tree.

The most efficient spanning tree can be found by a minimum spanning tree algorithm. The MST algorithm constructs a tree including every vertex, where the sum of the weights of all the edges in the tree is minimized. Prim’s algorithm [33] is used to construct the spanning tree by computing an edge with the least weight and adding it to the growing spanning tree. Prim’s algorithm is selected for our work since it is one of the fastest algorithms in dense graphs [34].

III-B4 Coverage Path Planning Algorithm

For the purpose of path planning, each block is logically divided into four parts, as depicted in Fig. 3, and the center of each part serves as the connection point along the constructed route. The construction of the coverage path depends on the number of adjacent blocks in a given direction, which can be categorized into three cases: (1) no adjacent block, (2) one adjacent block, or (3) multiple adjacent blocks. In the first case, the centers of the two parts are simply connected. In the second case, if the adjacent block is already connected to the current block, it is ignored. Otherwise, the two parts of the current block are connected to the adjacent parts of the neighboring block in the same direction. In the third case, the adjacent blocks are sorted in the same direction, and for each adjacent pair of blocks, a joint point is generated as the intersection point between two lines. The first line connects the midpoint between the two centers of the adjacent blocks to the center of the current block, and the second line connects the center of the two adjacent parts of the current block. Fig. 4 (d) illustrates these lines, and the relevant formulas are located in lines 16-20 of Algorithm 2. Finally, the adjacent parts of the adjacent blocks are connected to the generated joint points. Figs. 3 and 4 visually demonstrate the linking of blocks and the generation of the UGV route.

Denote $D_{d,b}$ as the list of neighbor blocks of the block $b$ in the direction $d$ , $Q_{c,b}$ as the coordinate of the center of the part $c\in\{TL,TR,BL,BR\}$ of the block $b$ , $P_{C,b}$ as the coordinate of the center of the block $b$ , $P_{L,b}$ as the coordinate of the left edge of the block $b$ , $N_{b}$ as the list of neighbor blocks of the block $b$ . The specific implementation is given in Algorithms 1 and 2.

Algorithm 2 Path Planning Algorithm

1:Inputs: List of blocks

B_{k}

, spanning tree

ST=\{E_{B_{i},B_{j}}|i\neq j,i,j=1,\dots,N_{B}\}

, starting point

p_{start}

2:Outputs: List of path elements

P

3:for

b\in B_{k}

b_{LEFT}\leftarrow D_{LEFT,b}

5: if

|b_{LEFT}|>0

then

6: if

|b_{LEFT}|=1

then

7: if

|D_{RIGHT,b^{1}_{LEFT}}|=1

then

8: Add path

(Q_{TL,b},Q_{TR,b^{1}_{LEFT}})

P

9: Add path

(Q_{BL,b},Q_{BR,b^{1}_{LEFT}})

P

10: end if

11: else

12: Sort

b_{LEFT}

x

in ascending order

13:

n\leftarrow|b_{LEFT}|

14: Add path

(Q_{TL,b},Q_{TR,b^{1}_{LEFT}})

P

15: for

i=1,\dots,n-1

16:

middle\leftarrow\frac{P_{C,b^{i}_{LEFT}}+P_{C,b^{i+1}_{LEFT}}}{2}

17:

center\leftarrow P_{C,b}

18:

joint_{x}\leftarrow P_{L,b}+\frac{BS_{b}}{4}

19:

joint_{y}\leftarrow middle_{y}+\frac{(joint_{x},middle_{x})*(center_{y}-middle_{y})}{center_{x}-middle_{y}}

20:

joint\leftarrow(joint_{x},joint_{y})

21: Add path

(Q_{BR,b^{i}_{LEFT}},joint)

P

22: Add path

(joint,Q_{TR,b^{i+1}_{LEFT}})

P

23: end for

24: Add path

(Q_{BL,b},Q_{BR,b^{n}_{LEFT}})

P

25: end if

26: else

27: Add path

(Q_{TL,b},Q_{BL,b})

P

28: end if

29: Do the same for the RIGHT, TOP, BOTTOM directions

30: Sort

P

in order of proximity to starting point

p_{start}

31:end for

III-B5 Predicting Time to Follow the Path

The turnaround time $\hat{T}$ and total angle difference $\hat{\theta}$ are estimated by:

\begin{split}&\hat{\theta}=\sum_{i=1}^{|P|}{|\alpha_{p_{i}}-\alpha_{p_{i+1}}|}\\ &\hat{T}=\frac{L}{v}+\frac{n_{t}*\hat{\theta}}{\omega},\\ \end{split}

(7)

where $\alpha_{p_{i}}$ is the orientation of the $i^{th}$ line $p$ with respect to the x-axis. The component $\frac{n_{t}*\hat{\theta}}{\omega}$ in the $\hat{T}$ equation is the leader’s rotational time. The total path length $\hat{L}$ and the coverage percentage $\hat{CP}$ are defined as in (3). At sharp corners, the leader only rotates around its heading; whereas, the followers modify their positions. The formation’s rotational motion at the center of the formation, therefore, can be approximated as that at the leader position.

III-B6 Optimising the Grid Cell Size to Meet the Time Budget

We use a binary search strategy with low memory and low complexity [35] to identify the appropriate grid cell size that produces a path length that can be followed within a given time budget. Binary search repeatedly divides the search interval in half of the lower $\underline{m}$ and the upper $\overline{m}$ bounds. If all estimated values are equal to the given metrics, the search is completed. The searching also stops when the lower $\underline{m}$ is greater than or equal to the upper $\overline{m}$ . Otherwise, if all estimated values of the current CS are more than the given metrics, the search narrows the interval in the upper half. Otherwise, the search recurs in the lower half. However, the maximum coverage percentage metric cannot be calculated accurately due to intermittent scanning values.

According to Algorithm 3, the binary search begins to compute $\hat{L}$ and $\hat{T}$ for every chosen grid cell size. It repeatedly checks until the two target values, which are less than $\hat{L}$ and $\hat{T}$ , are obtained. Then, the linear search resumes estimating $\hat{CP}$ for larger subsequent cell sizes until the final objective (maximum coverage percentage) is attained within 10 next cell size values by the linear search optimisation algorithm.

Algorithm 3 Grid Cell Size Optimisation Algorithm

1:Algorithm Parameters:

\hat{L}

\hat{T}

, and

\hat{CP}

2:Inputs: CSs,

\delta

3:Outputs:

\bar{CS}

4:Consider a cell size list with

n_{c}

element,

CS_{0},...,CS_{n_{c}-1}

, and three target metrics

\delta

m

is the pointer’s position while

\underline{m}

and

\overline{m}

are the left and right positions of search area.

5:Initialize

\underline{m}\leftarrow 0

\overline{m}\leftarrow n_{c}-1

i\leftarrow 0

6:while

\underline{m}<\overline{m}

7: Compute

\hat{L}_{m}

and

\hat{T}_{m}

at the middle position

m=\frac{\underline{m}+\overline{m}}{2}

8: if

\hat{L}_{m}>L_{max}

and

\hat{T}_{m}>T_{max}

then

9: Set

\underline{m}\leftarrow m+1

10: else if

\hat{L}_{m}=L_{max}

and

\hat{T}_{m}=T_{max}

then

11: Set

\underline{m},\overline{m}\leftarrow m

12: break

13: else

14: Set

\overline{m}\leftarrow m-1

15: end if

16:end while

17:

m\leftarrow\min(\underline{m},\overline{m})

18:while

m<n_{c}

19: Compute

\hat{CP}_{m}

at the

m

position.

20: if

\hat{CP}_{m}>CP_{max}

then

21:

CP_{max}\leftarrow\hat{CP}_{m}

22:

i\leftarrow 0

23: else

24:

i\leftarrow i+1

25: end if

26: if

i=10

then

27: break

28: end if

29:

m\leftarrow m+1

30:end while

31:Return

m

and

CS_{m}

III-B7 Formation Control

Each of the $n_{q}$ UGVs is given an identification number subscript $q_{0}$ , $q_{1}$ , etc. Each real UGV is assumed to be connected by a virtual spring system to its virtual leader, representing the physical interconnection between the UGVs and wireless communication. Referring to Fig. 5, orange lines represent the desired formation. Red circles indicate the follower UGVs, while green circles illustrate the virtual leader. Black circles are waypoints. Dashed green lines show the planned path. The black line is the spanning tree. Multiple UGVs ( $q_{1},...,and~{}q_{5}$ ) will be controlled to move in a formation whose pattern depends on the various block sizes. $q_{0}$ is treated as the virtual leader; the remaining UGVs are followers. Once virtual UGV $q_{0}$ and UGV $q_{i}$ ( $i\geq 1$ ) are linked together through virtual springs (VSs), spring forces are generated between them.

Based on the desired natural length vector $l_{o}$ and actual length vector $l_{a}$ of each VS, the common control rules of the VS method can be set.

The degree of difficulty in assigning a specific character for a UGV in the formation is defined as a character cost. The cost function equation is expressed as:

\phi_{ij}=\omega_{cx}|l_{a_{0ix}}|+\omega_{cy}|l_{a_{0iy}}|+\omega_{c\theta}|l_{a_{0i\theta}}|,\forall i,j\geq 1

(8)

where $\phi_{ij}$ is the cost of the UGV $i$ to obtain the $j^{th}$ goal role. $l_{a_{0ix}}$ and $l_{a_{0iy}}$ are the relative positions from the actual UGV position to the virtual leader position along the $x$ and $y$ axes. $l_{a_{0i\theta}}$ is the angular difference between the initial direction and the target direction. $\omega_{cx}$ , $\omega_{cy}$ , and $\omega_{c\theta}$ are weight constants.

After all costs for each character in the formation are calculated, the UGV with the maximum cost will be assigned a corresponding role in $fo$ , as shown in Fig. 5(a) and the following equation:

fo_{j}=i~{},~{}if~{}\phi_{ij}>\phi_{hj},\forall i\neq h,~{}i~{}\&~{}h\leq n_{q},\forall fo_{k}>i,k\neq j.

(9)

A relative position-based formation control method is derived here to maintain the desired formation shape for networked UGVs. The natural lengths of the springs $l_{0}=\{l_{01},...,l_{0n_{q}}\}$ are set according to the geometrical requirements for the desired formations.

Depending on the block size (BS), there are three formation shapes: the first is V-shaped (V-formation), the second is the U-shaped U-formation, and the third is a queuing (line) formation (Q-formation). The V-formation or U-formation is selected when the BS is greater than one. For a block size of 1, the Q-formation is used. Using the BS information and the formation type, the desired relative positions for each follower with respect to the virtual leader are computed to generate an obstacle avoidance formation, as described in Algorithm 4 of the Supplementary Material section. The control input to each ground vehicle is the resultant of the force vector generated by the virtual spring pairs connected to the UGVs.

III-B8 Path Tracking Algorithm

After the leader UGV’s coverage trajectory is computed, an online path planner outputs a series of way points $wp={wp_{0},...,wp_{k},...,wp_{n_{w}}}$ around the trajectory for the UGV navigation, where $n_{w}$ denotes the number of way points. Further, each UGV in the formation knows the remaining UGVs’ position information.

Using the virtual leader strategy, the spanning-tree path is shared among UGVs after the calculation process is completed. Based on the virtual vehicle tracking error of follower 1 (the closest follower, named $fo_{1}\in fo$ ), the proposed virtual velocity for the virtual vehicle $v_{q_{0}}$ is:

v_{q_{0}}=v\bigg{(}1-\min\big{(}\max(\frac{|{l_{a_{01}}}|-\underline{l_{q_{0}}}}{\overline{l_{q_{0}}}-\underline{l_{q_{0}}}},0),1\big{)}\bigg{)},

(10)

where $|{l_{a_{01}}}|$ represents the distance between the virtual leader and follower 1, $\overline{l_{q_{0}}},\underline{l_{q_{0}}}$ represent the maximum and minimum distance between the virtual leader and follower 1 to reach the minimum and maximum speed, respectively.

After every time step, the next way point $wp_{k+1}$ is produced:

wp_{k+1}:=wp_{k}+v_{q_{0}}~{}dt,

(11)

where $dt$ is the sample time.

Based on the Euler distance between the leader and the goal, the goal following force vector $F_{g}$ acting on the leader is derived as:

F_{g}=\omega_{g}*(wp_{k+1}-q_{p_{0}}),

(12)

where $\omega_{g}$ is the attractive force weight.

The proposed formation control approach is distributed to each UGV to guarantee that if any UGV fails the rest can adopt new roles and continue to track the virtual leader.

III-B9 Closest-Safe-Angle-Based Obstacle Avoidance

In order to prevent further collisions, a LiDAR sensor system is equipped on each UGV to measure distances ( $d_{o}$ ) and angles ( $\alpha$ ) from any surfaces. It works by emitting pulsed light waves in all directions and measuring how long it takes for them to bounce back off surrounding objects. To be convenient for further computation, the LiDAR information chain is converted to Cartesian coordinates as follows:

\alpha_{\iota}=\iota~{}\phi+\theta,\forall~{}0\degree\leq\alpha_{\iota}\leq 360\degree,

(13)

p^{\iota}_{o}=(x^{\iota}_{o},y^{\iota}_{o})=(d^{\iota}_{o}~{}\cos{\alpha_{\iota}},d^{\iota}_{o}~{}\sin{\alpha_{\iota}}),

(14)

where $\phi$ denotes an angular distance between measurements while $\theta$ illustrates the UGV heading and $\iota$ stands for the measurement step. $p^{\iota}_{o}=(x^{\iota}_{o},y^{\iota}_{o})$ is the obstacle position at the $\iota^{th}$ scan time. Additionally, $d^{\iota}_{o}$ is the $\iota^{th}$ relative distance measured at the $\iota^{th}$ $\alpha$ scanning angle.

When any of the UGVs or obstacles move or violate the avoidance radius $R_{a}v$ , the UGV’s angular width causing possible collisions called the blocked angles $\alpha_{b}$ , can be estimated and incorporated into the planner. To guarantee safe passage, the UGV is steered towards angles that are not blocked, called safe angles $\alpha_{s}$ .

A discrete list of all possible heading angles between 0 and 2 $\pi$ is generated based on the angle increment $\phi$ . For example, if $\phi=\frac{\pi}{180}$ , the list’s size $n$ will be 360. The algorithm then searches for all safe angles in the list: $0\leq\alpha^{\iota}_{s}<2\pi,\alpha^{\iota}_{s}\in\alpha\land 1\leq\iota\leq n$ .

We define a collection of obstacle cells that need to be avoided: $P^{\prime}={P^{k}_{o}:1\leq k\leq m}$ , the UGV circle’s radius $R$ , the UGV centre $O$ , the obstacle circle’s radius $R_{p}$ , the $P^{k}_{o}$ centre $C_{k}$ , and the relative angle $\eta^{k}$ between the centre line $OC_{k}$ and the $X-Axis$ . The two tangent lines between the circle of the obstacle $P^{k}_{o}$ and the UGV circle are constructed as shown in Fig. 6. Next, a set of two symmetric $k^{th}$ blocked angles [ $-\beta^{k}$ , $\beta^{k}$ ] with respect to the centre line $OC_{k}$ can be expressed as:

[-\beta^{k};\beta^{k}]=[-\arcsin{\frac{R+R_{p}}{OC_{k}}},\arcsin{\frac{R+R_{p}}{OC_{k}}}].

(15)

All global heading angles $\alpha$ located in the $\beta$ angle’s width are treated as blocked angles for the relevant UGV to transverse. The collection of angles blocked by the obstacle $P^{k}_{o}$ , named $\Delta_{k}$ , can be defined as:

\Delta_{k}:\{\alpha_{\iota}:\leq\alpha_{\iota}\leq\eta^{k}+\beta^{k},\alpha_{\iota}\in\alpha\}.

(16)

A complete list of the blocked angles $\Delta$ is obtained from:

\Delta=\cup^{k=1}_{m}\Delta_{k}.

(17)

Now, a list of safe angles $S$ can be established by excluding the list $D$ from the list $\alpha$ :

S=\alpha-\Delta.

(18)

If $S=\emptyset$ , the UGV would halt ( $V=0$ ) until any safe angle is scanned. Otherwise, the safe angle closest to the target angle $\alpha_{t}$ is selected as follows:

\alpha_{r}=min(|\alpha_{\iota}-\alpha_{t}|),\alpha_{\iota}\in S.

(19)

The outcome of our obstacle avoidance strategy is the UGV’s desired minimum orientation $\alpha_{r}$ . If no obstacles are predicted, the virtual force vector $(v_{x},v_{y})$ is fused from all force components (target force and spring force). Otherwise, this vector is turned towards the angle $\alpha_{r}$ in order to generate a new avoiding vector $v_{av}$ as in (20).

v_{av}=\begin{bmatrix}v_{avx}\\ v_{avy}\end{bmatrix}^{T}=\begin{bmatrix}v_{x}cos(\alpha_{r})-v_{y}sin(\alpha_{r})\\ v_{x}sin(\alpha_{r})+v_{y}cos(\alpha_{r})\end{bmatrix}^{T},

(20)

Based on (20), the control input of each UGV is computed as:

u=\begin{bmatrix}v\\ \omega\end{bmatrix}=\begin{bmatrix}v=|v_{av}|,\\ \omega=atan.(v_{av}),\end{bmatrix}

(21)

III-C Complexity Analysis

This section discusses the time complexity of the algorithm components mentioned above. Dynamic formation implementations demonstrate a worst-case time complexity of $O(n_{q}^{2})$ where $n_{q}$ stands for the number of mobile UGVs. This is because each physical UGV must inspect its position relative to the virtual leader and then transfer the relative position information to every other UGV to determine its role in the formation. The other significant element of our approach is the exchange of coverage matrices between UGVs. On receipt of such data from another UGV, each UGV must compare $C_{W}\times C_{H}$ cells to update its coverage and obstacle matrices. This process delivers time complexity $O(C_{W}C_{H})$ . Hence, the worst-case time complexity would be $O(n_{q}^{2})+O(C_{W}C_{H})$ if all agents exchange their coverage matrices and their own relative position information with all other agents. The algorithm is theoretically scalable to large formations and environments due to no exponential terms. In our case where there is a small number of UGVs, $O(C_{W}C_{H})$ is the dominant term.

IV Experiments

This section presents a range of experiments designed to evaluate the performance of our planning and prediction engine in various scenarios. In Section IV.A, we demonstrate the engine’s accuracy in predicting execution time by conducting an experiment with a simulated UGV. This experiment shows that the engine can efficiently find a suitable coverage path plan that meets a given time budget in only a short time.

In Section IV.B, we conduct a series of comparative experiments using our novel coverage path planning with formation control (CPPF), the coverage path planning with swarming, and the frontier-led swarming [30] on simulated Jackal UGVs. These experiments investigate the impact of UGV speed and time budget on path following performance while maintaining the group and order of the UGV formation. Section IV.C compares the multi-robot coverage path planning performance of CPPF and another method that uses a quadtree data structure without physical limits [36]. This comparison provides insights into the strengths and weaknesses of each approach to the optimal coverage path planning problem in cluttered environments. In the simulated experiments, we used the same settings (number of agents, maximum exploration time, path length, environmental setups, and initial locations) to guarantee a fair comparison between algorithms.

In Section IV.D, we describe experiments conducted on real Jackal UGVs in outdoor environments. These experiments aim to evaluate the engine’s performance in real-world conditions. Finally, we summarize the experiments in Section IV.E, providing an overview of the findings and their implications.

IV-A Experiment 1: Characterising Prediction Performance

IV-A1 Experiment 1 Setup

To verify the effectiveness of the prediction engine, a simulated, cluttered map of $25m\times 25m$ is used (see Fig. 1). Four static and complex-shaped obstacles (e.g., star, U-shaped, cylinder, and castle-shaped) were placed randomly. Cell sizes ranging from 0.75m to 3m are used to generate coverage paths and time-to-follow predictions. A single simulated UGV was then permitted to follow the path and the time taken compared to the prediction.

A single small Turtlebot UGV was used in these experiments so that we could examine the approximate path following performance without the complexity of formation control. The UGV was equipped with a 360-degree LiDAR sensor for reactive obstacle avoidance while path following. The forward speed $v$ of the UGV was initialised to 0.14 $m/s$ . The maximum turn rate when maneuvering was set to $0.7rad/s$ . Other parameters of this experiment are summarised in Table II. The i7-1260P CPU, Ubuntu Focal (20.04), ROS Noetic framework, and Python 3.8 are used to program and implement the planning approaches. The dynamic behavior of all vehicles was simulated in Gazebo. The sampling time of the whole system is set at 30Hz.

TABLE II: Experiment 1 parameters and their experimental values

Parameter	Description	Value
$\delta_{d}$	Distance tolerance	0.05m
$\delta_{\theta}$	Angle tolerance	15°
$\omega_{g}$	Target force weight	1.1

IV-A2 Experiment 1 Performance Metrics

Three criteria are chosen to evaluate the prediction performance of the proposed algorithms. They are:

•

PLD: Difference between the predicted and actual coverage path length,
•

TTD: Difference between the predicted and actual turnaround time (where turnaround time is defined as the time to complete following the path),
•

CPD: Difference between the predicted and actual coverage percentage.

IV-A3 Experiment 1 Prediction results

Tables III-IV show the results from the prediction module versus the Gazebo path following simulation. PLD was under 25m and TTD under 300s. Our observation of the UGV’s movement revealed that the need to slow down to perform turns was the main factor causing the error between predicted and actual performance.

Fig. 1 shows two examples of the predicted paths for different grid cell sizes. We can see that the larger grid cell size causes a coarser boundary around obstacles (shown in purple around the blue obstacles). Fig. 7 shows that this has the effect of reducing the total area that will be covered by the UGVs when the grid cell size is very large. Thus, while CPD was 0%, the actual coverage percent reduces with increasing grid cell size. This has the effect of meeting a tighter budget.

We also observed that the lower the grid cell size, the greater the computation time to make a prediction. For example, the computation time for the 0.75m cell size is approximately 6s, while that for cell size above 1.5m does not go above 1.5s.

TABLE III: Experiment 1: Predicted and actual path length (PPL and APL) metrics for different grid cell sizes.

Cell Size (m)	PPL (m)	APL (m)	PLD (m)
0.75	1228.5	1223.87	4.63
1.0	890.0	882.37	7.63
1.25	712.5	689.05	23.45
1.5	501	490.33	10.67
1.75	378	369.98	8.03
2.0	224	215.46	8.54
2.25	252	245.77	6.23
2.5	220	214.75	5.25
2.75	110	106.34	3.66
3.0	78	75.47	2.53

TABLE IV: Experiment 1: Predicted and actual turnaround time (PTT and ATT) metrics for different grid cell sizes.

Cell Size (m)	PTT (s)	ATT (s)	PLD (s)
0.75	12114.06	12161.37	-47.3
1.0	8640.02	8339.27	300.75
1.25	6440.17	6284.14	156.03
1.5	4386.41	4284.1	102.31
1.75	3180.21	3183.306	-3.09
2.0	1914.16	1846.89	67.27
2.25	2037.86	2019.13	18.74
2.5	1777.88	1756.58	21.3
2.75	884.45	876.92	7.53
3.0	597.53	626.17	-28.63

IV-B Experiments 2-4: Characterising Path Following Performance

In this section, we study the efficacy of our formation control system (denoted CPPF) in following a planned coverage path. We examine the impact of the top speed of the UGVs and the time budget on performance. We include three comparative algorithms: the first (denoted CPPS and shown in Algorithm 5) substitutes the leader-follower flexible formation strategy with a rule-based swarm strategy. A swarm can adapt its shape reactively to wider or narrower parts of the environment but does not require a leader. The second comparative approach (denoted FS [30]) uses a reactive, frontier-led swarming strategy with no planning.

Algorithm 5 Swarming waypoint follower

1:Algorithm Parameters: alignment radius scale factor

\mathcal{K}_{a}

, cohesion radius scale factor

\mathcal{K}_{c}

, separation radius scale factor

\mathcal{K}_{s}

2:Inputs: cell size

CS

, block size

BS

, number of UGVs

n_{q}

3:Outputs: swarm force

\sigma

4:The alignment radius

R_{a}

, cohesion radius

R_{c}

, separation radius

R_{s}

are determined as follows:

\begin{split}R_{a}&=\mathcal{K}_{a}\frac{CS~{}BS}{n_{q}}\\ R_{c}&=\mathcal{K}_{c}\frac{CS~{}BS}{n_{q}}\\ R_{s}&=\mathcal{K}_{s}\frac{CS~{}BS}{n_{q}}\end{split}

(22)

5:The swarm force

\sigma

is computed according to [30] with

R_{a}

R_{c}

R_{s}

IV-B1 Experiments 2-4 Simulation Setup

The simulation environments in these experiments are implemented on the same PC running Ubuntu, ROS Noetic, and Gazebo’s Jackal models as in Experiment 1. The simulated environment is a replica of the University of New South Wales Canberra campus, shown in Fig. 10, where light blue objects represent the real static obstacles. A UGV team comprising 5 Jackal mobile UGVs is used.

Each experiment is repeated 5 times. Mean and 95 $\%$ confidence intervals are reported where appropriate. The $P$ -Value from the Student T-test is used to distinguish the performance of different algorithms. If the variance of means between two sets is less than the expected $P$ value of 0.05 (i.e., 5 $\%$ ), we assume the results are statistically significant. The initial locations of the five vehicles are varied after every trial and given in Table V:

TABLE V: Experiments 2-4: Initial locations of UGVs

Trial Num	UGV 1	UGV 2	UGV 3
1	$[-23~{}0]\pm 0^{o}$	$[-26~{}-1]\pm 0^{o}$	$[-26~{}0]\pm 0^{o}$
2	$[-23.05~{}0.05]-10^{o}$	$[-26.05~{}-1.05]-10^{o}$	$[-26.05~{}0.05]-10^{o}$
3	$[-22.95~{}-0.05]+10^{o}$	$[-25.95~{}-0.95]+10^{o}$	$[-25.95~{}-0.05]+10^{o}$
4	$[-23.1~{}0.1]-20^{o}$	$[-26.1~{}-1.1]-20^{o}$	$[-26.1~{}1]-20^{o}$
5	$[-22.9~{}-0.1]+20^{o}$	$[-25.9~{}-0.9]+20^{o}$	$[-25.9~{}-0.1]+20^{o}$
Trial Num	UGV 4	UGV 5
1	$[-26~{}1]\pm 0^{o}$	$[-26~{}2]\pm 0^{o}$
2	$[-26.05~{}1.05]-10^{o}$	$[-26.05~{}2.05]-10^{o}$
3	$[-25.95~{}0.95]+10^{o}$	$[-25.95~{}1.95]+10^{o}$
4	$[-26.1~{}1.1]-20^{o}$	$[-26.1~{}2.1]-20^{o}$
5	$[-25.9~{}0.9]+20^{o}$	$[-25.9~{}1.9]+20^{o}$

Each UGV’s maximum linear speed is 2.0 $m/s$ , and the maximum angular velocity is set to $\pm 1.5rad/s$ . Table VI summarises all parameter settings for the simulation experiments.

Param	Unit	Description	Jackal
$\omega_{c}$	-	Weight for cohesion rule	0.27
$\omega_{a}$	-	Weight for alignment rule	1.05
$\omega_{s}$	-	Weight for separation rule	1.65
$\omega_{w}$	-	Wall Weight	1.1
$\omega_{g}$	-	Goal Weight	1.3
$\omega_{cx}$	-	X-Axis Weight	0.35
$\omega_{cy}$	-	Y-Axis Weight	0.35
$\omega_{c\theta}$	-	Yaw-Angle Weight	0.3
$R_{a}$	(m)	Alignment Radius	5
$R_{s}$	(m)	Separation Radius	1.1
$R_{c}$	(m)	Cohesion Radius	5
$R_{av}$	(m)	Avoiding Radius	0.3
$R_{d}$	(m)	Obstacle detection range	1.6
$k_{0}$	Spring coefficient	5.0
$\delta_{d}$	Distance tolerance	0.05m
$\delta_{\theta}$	Angle tolerance	15°
$\omega_{g}$	Target force weight	1.1
$v_{q_{0}}$	Virtual leader’s speed	1.1

TABLE VI: Experiments 2-4: Parameters of flexible formation control and the comparative frontier-led swarming algorithm.

The maximum coverage path length and coverage time ( $L_{max},T_{max}$ ) are 3500m and 30000s.

IV-B2 Experiments 2-4 Performance Metrics

In this series of experiments, the following metrics are examined:

•

Coverage percentage ( $CP$ ), the percentage of the environment visited by the UGVs
•

Path length ( $PL$ ), the length of the coverage path
•

Turnaround time ( $TT$ ), the time to follow the coverage path (or achieve 100% coverage in the case of the FS algorithm)
•

Coverage redundancy ( $CR$ ). Coverage redundancy (backtracking over areas already covered) can be calculated using (23) where $C_{repeated}$ is the average number of all repeated coverage cells.
•

Group ( $G$ ), how close together the UGVs are, defined as (24) where $N_{s}$ is the number of swarming UGVs, $N_{t}=5$ is the number of trials. $\bar{p}_{a}$ is the average position of the involved UGVs at the given time. We evaluate the group and order every 500 steps, as an average over the proceeding 150-time steps.
•

Order ( $O$ ), how well-aligned UGVs are in terms of both speed and direction as defined in 25) where $\bar{\varepsilon}_{a}$ is the average velocity of the involved UGVs at the given time. We also evaluate ordering every 500 steps as an average over the proceeding 150-time steps.

CR=\frac{|C_{repeated}|}{|C_{free}|}\times 100\%,

(23)

G=\frac{\sum_{i=0}^{N_{t}}\sum_{t=T_{0}}^{TT}\frac{\sum_{i=1}^{N_{s}}||p^{i}-\bar{p}_{a}||_{2}}{TT-t}}{N_{s}~{}N_{t}},

(24)

O=\frac{\sum_{i=0}^{N_{t}}\sum_{t=T_{0}}^{TT}\frac{\sum_{i=1}^{N_{s}}||\varepsilon^{i}-\bar{\varepsilon}_{a}||_{2}}{TT-t}}{N_{s}~{}N_{t}},

(25)

IV-B3 Experiment 2 Group and Order results

First, we examine the ability of all three algorithms to keep the UGVs together and heading in the same direction and at the same speed. Fig. 8 indicates that all the approaches maintain a reasonable level of grouping and ordering relative to the search space size. The FS approach maintains a tighter formation than the CPPF and CPPS approaches in all experimental settings. The difference in group metrics is statistically significant at the 95% confidence level. This is likely due to the influence of the block-size switches. However, grouping and order are still maintained by the CPPF and CPPS methods, and the group and order metrics do not change significantly over time. The order metrics are similar for all three approaches. This makes sense because the swarming and formation approaches should both encourage ordering.

IV-B4 Experiment 3 Coverage performances with three different speed caps

Three maximum vehicle speeds are set as 0.45m/s (low), 0.7m/s (medium), 0.95m/s (high) in this experiment. This permits us to examine whether high speeds cause the UGVs to miss covering parts of the environment. Based on the default coverage path length budget, turnaround time budget, and the maximum linear and rotation speeds, the optimisation algorithm proposes a grid cell size of 2.25m that should result in a coverage path with $\hat{L}$ of 2007.34m and $\hat{T}$ of 5039.04s for the low speed; 3445.92s for the medium speed and 2991.28s for the high speed. $\hat{CP}$ should be 100%. The LiDAR sensor’s observation range is 4m.

Fig. 9 shows the turnaround time, path length, coverage percentage, and coverage redundancy metrics.

As expected, the higher the maximum vehicle speed, the shorter the turnaround time is. Using CPPF and FS, the TTD is below 300s, while TTD for CPPS is approximately double. CPPF and FS strategies were also able to cover the whole region ( $78.19\pm 0.00\%$ direct $CP$ plus $21.81\pm 0.00\%$ indirect $CP$ ), whereas the CPPS achieved an incomplete coverage percent ( $77.37\pm 0.32\%$ for the direct coverage and $21.99\pm 0.00\%$ for the indirect coverage). This is because CPPS uses an organic formation control strategy that is unable to backtrack.

CPPF has a PLD of under 90m, while FS has the largest PLD. This is because FS (which does not include a planning component) often needs to backtrack to cover a missed area. The differences in CR for the comparative algorithms are insignificant at the 95% confidence level, except in the high-speed case, where the FS produces the highest redundant coverage percentage, namely, $2.1\pm 0.06\%$ . While the path planning-based coverage methods generate the non-backtracking paths, the frontier search technique guides the physical UGVs to track frontier cells, leading to an increase in back-tracking the covered areas.

IV-B5 Experiment 4 Performance with different time and path length budgets

All algorithms were also run with three different time and path length budgets: (1) high budgets of $[4000s,2500m]$ (more than enough time for all algorithms to achieve 100% coverage), (2) tight budgets of $[2800s,1800m]$ (barely enough time for one of the algorithms to finish), (3) very tight budgets of $[1200s,800m]$ (no algorithm will be able to achieve optimal solution-approximately 77.15% map covered).

The planner produces the paths shown in Fig. 10 for each of these cases. In these figures, blue polygons show the real obstacle positions. Purple cells indicate regions denoted as obstacles as a result of different grid cell size recommendations. Yellow cells represent open areas. All UGVs’ maximum speed in all experiments is set at 0.7m/s.

In general, Fig. 11 demonstrates that the CPPF strategy presents minimal offsets from the desired coverage parameters (such as the maximum coverage percentage, total path length, and coverage time) in both cases. The FS and CPPF exhibit quite similar metrics for the very tight and high budgets scenarios. In contrast, the FS has a superior $TT$ figure for the tight-budget case, with the lowest $PL$ , and $CR$ values. This difference is statistically significant at the 95% confidence level. Moreover, there are no significant differences in the $TT$ and $CR$ parameters when high budgets are allowed.

As shown in Fig. 11(c), it is interesting to see that both strategies only achieve partial (direct) coverage in the tight and very tight-budget scenarios because all obstacle boundaries are not observed by the LiDAR sensor (its observation range is smaller than the grid cell size). They produce the same direct $CPs$ of only $45\pm 0.00\%$ and $68.46\pm 0.00\%$ . However, they generate full coverage, including direct and indirect coverage, when a high budget is employed. The $CP$ metric of the CFFS is $0.2\%$ lower than those of the CPPF and FS. The reason is that the higher the desired budget is set, the smaller the grid cell size is. In the high-budget case, the grid cell size of 2.25m is significantly less than LiDAR’s observation range of 4m.

IV-C Experiment 5: Comparison with a Quadtree Data Structure

The spanning tree produced by the multi-resolution block structure used in CPPF can have a significant impact on coverage results. To assess this impact, we compared the performance of MCPP in CPPF with that of a conventional homogeneous quadtree algorithm, where all blocks are of uniform size (0.45m), using a 25m $\times$ 25m grid map containing oddly shaped obstacles as illustrated in Fig. 1. Complex-shaped obstacles can cause a block to have multiple neighbor blocks in the same direction, resulting in incomplete coverage routes (green line) and discontinuous minimum spanning trees (black line) in the quadtree method. However, our MCPP algorithm resolves this issue in the third case of Algorithm 2 (Lines 12-24). In this case, when obstacle cells exist between two adjacent blocks, the joint point still satisfies the condition of being within the intersection region between the considered block and the triangle formed by the three centers of the considered block and the two adjacent blocks (see Fig. 12).

As observed in Figs. 13 and 14, blocks B2 and B3 are located on top of block B1, as determined by the MST. The quadtree algorithm is then applied to link the two upper parts of B1 with the two lower parts of B2, thereby precluding the possibility of establishing an additional connection with B3. The disruption of the connection between B1 and B3 would result in the loss of the remaining branch of the MST. In the absence of obstacles on the map, the minimum spanning tree will exhibit, at most, a single connection between any two blocks in any given direction, and the seamless continuity of the path is ensured by the quadtree algorithm. Our proposed algorithm establishes a path that connects B1 to both B2 and B3, resulting in a path that guarantees that all parts of the MST are followed to form the robot path.

Furthermore, although all estimates, such as maximum coverage percentage, total path length, and coverage time, can be well achieved by both methods when combined with our formation control solution, as demonstrated in Experiments 1-4, the coverage percentage obtained by running the quadtree algorithm is lower than our spanning tree method with modified rules. This is because the multi-robot team cannot move through the blocks, such as the eleven yellow blocks shown in Fig. 12, to update the coverage information due to the absence of the coverage paths.

IV-D Experiment 6: Outdoor experiments using Jackals and DGPS positioning

IV-D1 Experiment 6 Outdoor Setup

To verify the effectiveness of the algorithms on real UGVs, a $15.25m\times 24.4m$ outdoor setup was used involving several arbitrary-shaped obstacles (see Fig. 15) and a team of three Jackal UGVs. There were three significant obstacles: a hut, a shipping container, and a piece of equipment covered with a metal mesh cage. Several steel posts were on the terrain, each with a Y cross-section. The diameter of a post (assuming a circular cross-section) was roughly 6cm. In addition to these fixed obstacles, a no-go zone for the UGVs was included to prevent the UGVs from crossing over known rabbit burrows.

The Jackal UGVs used in our experiments are fitted with a differential GPS (DGPS) rover module, Inertial Measurement Unit (IMU) and a SICK LMS-111 LiDAR. The LiDAR’s observation range is 4m. Incoming GPS rover signals and the internal IMU sensor were read at a rate of 10Hz. The DGPS base station was stationary and transmits DGPS corrections to the rovers. The sampling time of the whole system was set at 30Hz. The role of the ROS Master is to enable individual ROS nodes to locate one another. Unlike the traditional implementations that run the ROS Master on only a single base station, the ROS Master was implemented on each UGV to facilitate ROS communications and improve the system’s robustness.

Five trials were conducted. The problem constraints given were the maximum path length $L_{max}$ of 100m, and maximum coverage time $T_{max}$ of 600s. For each test, the UGVs were initially setup in a V-formation close to the bottom-left corner at GPS-measured coordinates of [9.53, 15.58, 0]m, with the formation facing down towards the bottom of the map. Some slight variations in initial UGV locations were introduced at the start of each test. The initial locations of the three vehicles are given in the columns of (26). The virtual leader’s forward speed $v$ is initialised to 0.2 $m/s$ . The maximum turn rate when maneuvering is set to $0.73rad/s$ . Other parameters of our experiments are summarised in Table VII.

q_{p_{i}}(0)=\begin{bmatrix}9.23&10.01&8.07\\ 15.69&18.4&18.32\\ 0&0&0\end{bmatrix}

(26)

TABLE VII: Experiment 6: Parameters of leader-follower formation control

Parameter	Description	Jackal
$k_{0}$	Spring coefficient	3.9
$\zeta$	Vision range of LiDAR sensor	4.0m
$\delta_{d}$	Distance tolerance	0.3m
$\delta_{\theta}$	Angle tolerance	15°
$R_{av}$	Avoidance radius	0.32m
$\omega_{g}$	Target force weight	1.1

IV-D2 Experiment 6 Performance metrics

The metrics used in this experiment were: (1) difference between the actual and predicted path length (PLD); (2) difference between the actual and predicted turnaround time (TTD); (3) Coverage percentage (CP); (4) coverage redundancy (CR); (5) group (G); and (6) order (O).

IV-D3 Experiment 6 Results

The experimental test performed can be viewed in the following videos: https://youtu.be/z4kK6OnXXg8. The prediction engine took 2s to produce the recommended path shown in Fig. 16. It recommends a grid cell size of $3.05m$ , and predicts the path length $\hat{L}$ of 96.86m, coverage time $\hat{T}$ of 575.04s, and coverage percentage $CP$ of 100 $\%$ .

As shown in the video, the map is entirely covered (namely $90\pm 0.00\%$ for the direct coverage and $10\pm 0.00\%$ for the indirect coverage) after $586.17\pm 1.2$ s (mean over five tests $\pm$ standard deviation). The travelled path length is approximately $103.03\pm 0.9$ m. These results reflect a PLD of 6.17m and TTD of 11.13s. These results are slightly higher than those from our simulations. We observed this is due to the variations in the initial UGV position among five tests and the actual tracking errors caused by the terrain’s uneven round with thick grass and the DGPS and IMU sensor errors. On the other hand, the advantage of using the spanning tree to design the complete coverage paths can be seen through the low CR of $7.33\pm 1.41$ when the back-tracking routes are eliminated, and only several cells at the corners are covered repeatedly.

As shown in Fig. 17, there are slight fluctuations in the actual paths. This is due to three factors. First, we measured GPS static position errors of 0.08m that contribute to some deviations in movement. Secondly, as the robots switch between V, U and queuing formations, there is some distortion in the paths. Finally, inter-UGV collision avoidance exerts influence on the robots causing deviations from the path.

However, under the control of our formation and role assignment strategies, three Jackal follower UGVs switch positions and effectively form the desired formations when tracking the virtual leader’s motion and avoiding obstacles along the designed coverage path. The UGV behaviors exposed in the real-time environment are similar to those obtained in simulations. As a result, the proposed methods yield reasonable $G$ and $O$ figures, although there are short periods where there is high variance in the grouping. This occurs when one UGV strays away from the others (e.g. as a result of a brief increase in GPS error). The system self-corrects when GPS is re-acquired. Additionally, there are no significant variations of these two variables over time (see Fig. 18).

Besides stationary obstacles, dynamic obstacles such as humans can move unpredictably and obstruct a flock’s movement along a planned path. However, by combining the MCPP method with the closest safe angle-based obstacle avoidance, our robot team was able to successfully navigate through a challenging environment without any communication. In a video demonstration (available at https://youtu.be/UjihyOj-VCU), the relevant UGVs quickly adjusted their path when the human obstacle violated their obstacle avoidance radius, safely steering towards the nearest safe area and then resuming their intended path. Furthermore, during formation switches, the UGVs were able to avoid mutual collisions. This work addresses the limitations of existing coverage path planning methods by incorporating obstacle avoidance techniques, which enable robots to navigate safely in a dynamic environment [37].

IV-E Summary of Experiments

To validate the effectiveness and efficiency of our algorithm for multi-agent systems, we conducted a comprehensive set of experiments in both simulated and real-world environments, with five trials for each experiment. The obtained results showed that our algorithm yielded the best coverage accuracy rate without collisions and a coverage redundancy rate of only 2.7% when all physical limits were met. These results demonstrate our algorithm’s high accuracy and efficiency in ideal conditions.

We encountered some limitations in the real-time field tests due to hardware and environmental factors and mobile obstacles. However, even under these challenging conditions, we still achieved similar results to those obtained in the simulated experiments. This indicates that our algorithm is robust and can perform and adapt well to changing real-world scenarios.

The significance of these results lies in the potential applications of multi-agent systems, such as search and rescue operations or environmental monitoring. Our algorithm provides an efficient and reliable method for coordinating multiple agents to explore a dynamically changing area with minimal overlap and using minimal resources.

V Conclusion

This paper has described a novel approach to MCPP in a dynamic environment. The CPPF algorithm produces a coverage path that maximises the area of the environment that will be visited while meeting given time and path length budgets. Further, it permits the detection of unknown obstacles simultaneously with the steering of the mobile robot to avoid collisions. We demonstrated this algorithm on simulated and real Jackal UGVs. We provided a range of statistics showing the performance of the prediction engine and the path-following algorithms. We saw that: We saw that:

•

The planner is able to recommend suitable grid cell sizes to meet given time and path length budgets within approximately 10s.
•

Path length predictions are reasonably accurate for both simulated and real UGVs. The longer the necessary path, the lower the path length prediction error, with as little as 0.4% error on a 1km path.
•

Turnaround time predictions range in accuracy from a few seconds over 100 metres to up to 5 minutes discrepancy on a 1km path.
•

Formation-based path following achieves comparable coverage performance to a state-of-the-art reactive swarming approach but offers the guarantee of a time and path length prediction.
•

It is feasible to use the algorithm on real UGVs in an outdoor setting.

The possibilities for future work in this area are rich. The proposed method constructs the optimal coverage path for every UGV using the MST such that the union of all paths generates a full coverage of the terrain. However, these paths are static, and the coverage is performed in static environments where the obstacles do not move. In future work, we will use a rapidly-exploring random tree (RRT) algorithm for local path re-planning to update the current path to achieve mobile obstacle avoidance in dynamic environments.

VI Acknowledgement

This work was supported by the Australian Defence Science and Technology Group (DSTG) under grant No. 9729.

References

[1] G. Huang, X. Yuan, K. Shi, Z. Liu, and X. Wu. A 3-d multi-object path planning method for electric vehicle considering the energy consumption and distance. IEEE Transactions on Intelligent Transportation Systems, 2021.
[2] E. Galceran and M. Carreras. A survey on coverage path planning for robotics. Robotics and Autonomous systems, 61(12):1258–1276, 2013.
[3] R. J. Wai and A. S. Prasetia. Adaptive neural network control and optimal path planning of UAV surveillance system with energy consumption prediction. IEEE Access, 7:126137–126153, 2019.
[4] H. X. Pham, H. M. La, D. Feil-Seifer, and M. Deans. A distributed control framework for a team of unmanned aerial vehicles for dynamic wildfire tracking. In 2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pages 6648–6653. IEEE, 2017.
[5] Y. Wang, J. Xu, Q. Liu, Y. Zhang, and J. Yang. Path planning of seeding robot based on improved ant colony algorithm. In Proceedings of 2021 Chinese Intelligent Automation Conference, pages 31–37. Springer, 2022.
[6] C. Liu, S. Zhang, and A. Akbar. Ground feature oriented path planning for unmanned aerial vehicle mapping. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 12(4):1175–1187, 2019.
[7] F. Niroui, K. Zhang, Z. Kashino, and G. Nejat. Deep reinforcement learning robot for search and rescue applications: Exploration in unknown cluttered environments. IEEE Robotics and Automation Letters, 4(2):610–617, 2019.
[8] R. Almadhoun, T. Taha, L. Seneviratne, and Y. Zweiri. A survey on multi-robot coverage path planning for model reconstruction and mapping. SN Applied Sciences, 1(8):1–24, 2019.
[9] L. J. Colombo and H. G. de Marina. Forced variational integrators for the formation control of multiagent systems. IEEE Transactions on Control of Network Systems, 8(3):1336–1347, 2021.
[10] V. P. Tran, M. A. Mabrok, M. A. Garratt, and I. R. Petersen. Hybrid adaptive negative imaginary-neural-fuzzy control with model identification for a quadrotor. IFAC Journal of Systems and Control, 16:100156, 2021.
[11] V. P. Tran, M. A. Garratt, and I. R. Petersen. Multi-vehicle formation control and obstacle avoidance using negative-imaginary systems theory. IFAC Journal of Systems and Control, 15:100117, 2021.
[12] J. M. Palacios-Gasós, E. Montijano, C. Sagüés, and S. Llorente. Distributed coverage estimation and control for multirobot persistent tasks. IEEE transactions on Robotics, 32(6):1444–1460, 2016.
[13] J. M. Palacios-Gasós, D. Tardioli, E. Montijano, and C. Sagüés. Equitable persistent coverage of non-convex environments with graph-based planning. The International Journal of Robotics Research, 38(14):1674–1694, 2019.
[14] S. Manjanna, A. Q. Li, R. N. Smith, I. Rekleitis, and G. Dudek. Heterogeneous multi-robot system for exploration and strategic water sampling. In 2018 IEEE International Conference on Robotics and Automation (ICRA), pages 1–8. IEEE, 2018.
[15] V. P. Tran, M. A. Garratt, K. Kasmarik, and S. G. Anavatti. Dynamic frontier-led swarming: Multi-robot repeated coverage in dynamic environments. IEEE/CAA Journal of Automatica Sinica, 10(3):1–16, 2023.
[16] W. Dong, S. Liu, Y. Ding, X. Sheng, and X. Zhu. An artificially weighted spanning tree coverage algorithm for decentralized flying robots. IEEE Transactions on Automation Science and Engineering, 17(4):1689–1698, 2020.
[17] X. X. Shao, Y. J. Gong, Z. H. Zhan, and J. Zhang. Bipartite cooperative coevolution for energy-aware coverage path planning of uavs. IEEE Transactions on Artificial Intelligence, 3(1):29–42, 2021.
[18] J. Chen, C. Du, Y. Zhang, P. Han, and W. Wei. A clustering-based coverage path planning method for autonomous heterogeneous uavs. IEEE Transactions on Intelligent Transportation Systems, 2021.
[19] C. Gao, Y. Kou, Z. Li, A. Xu, Y. Li, and Y. Chang. Optimal multirobot coverage path planning: ideal-shaped spanning tree. Mathematical Problems in Engineering, 2018, 2018.
[20] G. E. Jan, C. Luo, H. T. Lin, and K. Fung. Complete area coverage path-planning with arbitrary shape obstacles. Journal of Automation and Control Engineering Vol, 7(2), 2019.
[21] L. Zhou, V. Tzoumas, G. J. Pappas, and P. Tokekar. Distributed attack-robust submodular maximization for multirobot planning. IEEE Transactions on Robotics, 2022.
[22] V. P. Tran, M. Garratt, and I. R. Petersen. Switching time-invariant formation control of a collaborative multi-agent system using negative imaginary systems theory. Control Engineering Practice, 95:104245, 2020.
[23] A. Gorbenko and V. Popov. The multi-robot forest coverage for weighted terrain1. Journal of Ambient Intelligence and Smart Environments, 7(6):835–847, 2015.
[24] G. Q. Gao and B. Xin. A-stc: Auction-based spanning tree coverage algorithm formotion planning of cooperative robots. Frontiers of Information Technology & Electronic Engineering, 20(1):18–31, 2019.
[25] Y. Ma, Y. Zhao, Z. Li, H. Bi, J. Wang, R. Malekian, and M. A. Sotelo. CCIBA*: An improved BA* based collaborative coverage path planning method for multiple unmanned surface mapping vehicles. IEEE Transactions on Intelligent Transportation Systems, 2022.
[26] K. O. Ellefsen, H. A. Lepikson, and J. C. Albiez. Multi-objective coverage path planning: Enabling automated inspection of complex, real-world structures. Applied Soft Computing, 61:264–282, 2017.
[27] K. O. Ellefsen, H. A. Lepikson, and J. C. Albiez. Planning inspection paths through evolutionary multi-objective optimization. In Proceedings of the Genetic and Evolutionary Computation Conference 2016, pages 893–900, 2016.
[28] B. Zhou, J. Pan, F. Gao, and S. Shen. Raptor: Robust and perception-aware trajectory replanning for quadrotor fast flight. IEEE Transactions on Robotics, 37(6):1992–2009, 2021.
[29] S. Wen, Z. Wen, D. Zhang, H. Zhang, and T. Wang. A multi-robot path-planning algorithm for autonomous navigation using meta-reinforcement learning based on transfer learning. Applied Soft Computing, page 107605, 2021.
[30] V. P. Tran, M. A. Garratt, K. Kasmarik, S. G. Anavatti, and S. Abpeikar. Frontier-led swarming: Robust multi-robot coverage of unknown environments. Swarm and Evolutionary Computation, 75:101171, 2022.
[31] A. Aghaeeyan, F. Abdollahi, and H. A. Talebi. UAV–UGVs cooperation: With a moving center based trajectory. Robotics and Autonomous Systems, 63:1–9, 2015.
[32] J. H. Clark. A fast algorithm for rendering parametric surfaces. In Computer Graphics (SIGGRAPH’79 Proceedings), volume 13, pages 7–12. Citeseer, 1979.
[33] H. J. Greenberg. Greedy algorithms for minimum spanning tree. University of Colorado at Denver, 1998.
[34] F. Huang, P. Gao, and Y. Wang. Comparison of prim and kruskal on shanghai and shenzhen 300 index hierarchical structure tree. In 2009 International Conference on Web Information Systems and Mining, pages 237–241. IEEE, 2009.
[35] K. K. Huang and D. Q. Dai. A new on-board image codec based on binary tree with adaptive scanning order in scan-based mode. IEEE Transactions on Geoscience and Remote Sensing, 50(10):3737–3750, 2012.
[36] X. Huang, M. Sun, H. Zhou, and S. Liu. A multi-robot coverage path planning algorithm for the environment with multiple land cover types. IEEE Access, 8:198101–198117, 2020.
[37] D. Jang, J. Yoo, C. Y. Son, D. Kim, and H. J. Kim. Multi-robot active sensing and environmental model learning with distributed gaussian process. IEEE Robotics and Automation Letters, 5(4):5905–5912, 2020.