Dynamic Elliptical Shaping Control for Swarm Robots

Solving the robotic swarm coverage problem for an elliptical area has various applications for exploring novel environments. Solutions for this problem should cover a specified ellipse and seamlessly adapt to changing numbers of robots. Previous solutions used techniques such as formation control, vector fields, and neural networks. While these techniques were successful, they all lacked one or more of the three key tenants of swarm elliptical attraction: complete coverage of an ellipse with commandable parameters, simplicity for scaling in the number of robots, and adaptive sizing. Additionally, no previous work presented guidelines for ensuring that the swarms could successfully and safely converge to the commanded ellipse without collisions. In contrast, this work presents a novel swarm elliptical attraction behavior with all three key tenants with guidelines for ellipse and swarm parameter selection. First, a new Lyapunov stable elliptical attraction behavior for Reactive Particle Swarms is presented. The behavior commands robots to cover the entire ellipse area for a specific semimajor axis, eccentricity, and orientation. Additionally, dynamic interagent spacing naturally ensures coverage for different numbers of robots. Second, the work presents a novel adaptive sizing algorithm that varies the ellipse’s semimajor axis based on the swarm state. The adaptive sizing algorithm specifies the eccentricity and orientation using time-varying functions. Third, guidelines for selecting the number of robots, commanded ellipse area, obstacle avoidance distance, and robot communication range that allow for successful aggregation to the commanded ellipse are presented. All three of the results are verified using simulation and hardware-in-the-loop trials.


I. INTRODUCTION
Swarm robotics uses the decentralized control of multiple agents to perform global behaviors through local individual decisions and are characterized by their simplicity, the interchangeability of agents, and the ability to scale in number [1]. There are a variety of subcategories of swarms, including leader-follower [2], consensus [3], [4], Particle Swarm Optimization [5], and Reactive Particle Swarms (RPS). This work uses RPS for its simplicity and reactivity to the environment and the swarm state. The RPS Control Architecture was developed in [6].
Multi-agent robotic systems naturally lend themselves to the distribution problem, which attempts to space robots throughout a specific area. Various swarm robotic solutions have been developed for these problems, from simple nearest neighbor dispersal [7] to artificial potential fields with formation control [8] to neural networks [9]. Additionally, intricate shapes beyond the traditional simple geometric shapes can be covered using bounding functions and superposition [10], [11].
An exciting subset of the distribution problem is the elliptical shaping problem. There are three key tenants to swarm elliptical attraction. First, the swarm should cover the entire elliptical area with specifiable ellipse parameters of semi-major axis, eccentricity, and orientation. This includes the perimeter of the ellipse and the ends when highly eccentric. Second, as with all swarm behaviors, the algorithm should be simple enough to scale to a high number of robots while also dynamically handling the addition or loss of robots. Third, the algorithm should be able to size the commanded ellipse adaptively, responding to the swarm state and the environment [8], [9]. Elliptical shaping can benefit the exploration of novel environments, including giving better sensor measurements across features of interest and adaptive sizing allows the swarm to scale and rotate to avoid obstacles [12].
Previous work demonstrated swarm attraction to a virtual elliptical area with set parameters and inter-agent positions set using explicit formation control [8], [13]. Using formation control in this manner is challenging because the formation needs to be recalculated every time the ellipse parameters or the number of robots in the swarm change. Other elliptical attraction behaviors use the superposition of vector fields to attract the robots and distribute them throughout the ellipse [14], [15]. While these methods have controllable ellipse parameters, they do not allow for adaptive sizing, and some use complex limit functions to force robots to distribute throughout the ellipse. Some work has shown varying the eccentricity of the commanded ellipse using artificial potential fields to navigate obstacles in a workspace [16]. This work has coverage around the circumference of the ellipse, not the interior region. Other work solves the estimation and formation control problem when communication ranges are very limited [17]. While this work solves the harder problem of estimation in addition to elliptical shaping, it has only been shown for a small number of robots and was unable to achieve good coverage of the ellipse area.
Reference [18] presents a discrete-time kinematic model for aggregation in a multi-dimensional space. Their algorithm bounds the swarm to elliptical regions in 2D and ellipsoids in 3D. Within the commandable bounding region, the swarm only preserves inter-agent connectivity and does not cover the region. Reference [9] uses a neural network to change the shape of the swarm to pass through narrow channels. However, the algorithm does not allow specific commanding of the ellipse parameters and does not have good coverage of the ends of the ellipse for highly eccentric shapes [9]. Voronoi Tessellations are another popular solution but can be slow to respond to changing ellipse parameters [19]. The slow convergence issue was improved using feedforward action but is still relatively computationally extensive for large numbers of robots [20].
Reference [12] is the most complete elliptical distribution algorithm in literature, allowing for scaling and rotations of a specified ellipse. The inter-robot distance is explicitly set and scaled with the ellipse to ensure decent coverage of the ellipse area for that number of robots. Because this distance is constant under scaling, it does not guarantee that the swarm will completely cover the ellipse, especially when robots are added or removed. Additionally, [12] does not provide guidelines for selecting the inter-robot distance.
Overall, previous methods can elliptically attract the swarm. However, none have all the desired elliptical attraction functions of commandable parameters, simplicity for scaling in the number of robots, and adaptive sizing. We chose to use the RPS architecture to develop an elliptical attraction function that meets the three tenants of elliptical attraction. This choice was motivated by the RPS architecture's natural handling of changes in the number of swarm agents and the emphasis on simple, scalable algorithms. Additionally, the RPS architecture is good at dynamic control, which should yield quick performance when tracking changing ellipse parameters. The proposed elliptical attraction function provides complete coverage of the ellipse area with specified parameters. The commanded ellipse can adaptively size, and the robot density is naturally managed to accommodate changing ellipse parameters and variations in the number of robots. Additionally, our elliptical attraction function ensures coverage of the entire ellipse, including the circumference.
The obstacle avoidance range must be preserved to avoid interagent collisions. However, it might be impossible to adhere to that obstacle avoidance range while fitting the swarm in the commanded ellipse. None of the previously mentioned works provide guidelines for ensuring that the swarm can fit in that desired ellipse which is vital for the widespread implementation of these algorithms. This paper has three major contributions. First, the work presents a novel Lyapunov stable RPS elliptical attraction behavior with coverage of the ellipse area and circumference and adheres well to the commanded ellipse parameters. Second, a new adaptive sizing algorithm for the elliptical attraction behavior is presented. The algorithm scales the semi-major axis as a function of the swarm state while ensuring successful elliptical aggregation. Third, novel guidelines for swarm and ellipse parameter selection are outlined to ensure successful elliptical attraction. These guidelines relate the ellipse area, communication range, obstacle avoidance range, and the number of robots for successful aggregation. All results are verified using simulated and hardware-in-theloop case studies.
The remainder of the paper is subdivided as follows. Section II summarizes the RPS architecture. Sections III and IV respectively present the RPS elliptical attraction behavior and the adaptive sizing algorithm. Section V outlines the design processes for the simulation and hardware case studies. Then Section VI outlines guidelines for swarm and ellipse parameter selection to ensure successful ellipse aggregation. Finally, Sections VII and VIII present a series of hardware verification trials for the elliptical attraction behavior and adaptive sizing algorithm.

II. RPS ARCHITECTURE SUMMARY
This paper uses the planar Reactive Particle Swarm (RPS) control architecture to design an elliptical attraction behavior. RPSs are homogenous swarms with scalable, lightweight, and reactive behaviors that do not require global knowledge, use consensus filtering, or have a hierarchy of robots. RPSs respond to the current swarm state with simple algorithms allowing for reactive control. RPS have broad applications, including various aggregation behaviors, collision avoidance, coverage, formation control, and gradient estimation [6], [21]. The RPS control architecture is a unifying control architecture that allows for the rapid development of novel swarm behaviors. The RPS architecture is summarized in [6], and Fig. 1 shows the architecture block diagram. In RPSs, there are two layers of control, the upper composite behavior layer and the lower base behavior layer. The composite behavior command velocity for robot i, ⃗ v ci , is a weighted sum of p base behavior command velocities, ⃗ v bi , and q external input velocities ⃗ u ui , as in (1). The scaling constants k are constrained as in (2).
The b th base behavior command velocity for robot i, ⃗ v bi , is calculated as in (3) using the pair ( calculates using a reference vector using the weighing function f b . Then, the pre-multiplication by K (θ b ) rotates the reference vector by θ b about the z-axis to calculate the behavior command velocity.
To calculate the reference vector, the relative position vectors pointing from robot i to all N robots in the swarm are calculated and vertically augmented together as ⃗ r i ∈ R 2N . Robots outside the communication range are assumed to be infinitely far away. The matrix D i normalizes the individual relative position vectors to relative unit vectors. Normalizing the vectors allows the behavior to be purely determined by the weighting function. Then C i (d l ) uses d l , the artificial communication range (ACR), to limit which robots robot i can use information from. This is done by zeroing the relative unit vectors for robots that are further than d l away from robot i. The ACR should be less than or equal to the actual communication range of the robots. Using an ACR smaller than the actual communication range allows the behavior to both avoid intermittent communication at the limit of the communication range and control the local region of information used.
Next, the weighted sum of the ACR limited relative unit vectors is calculated by multiplying by is a weighting function that returns for a pair of robots i and j a 2-by-2 weighting matrix. The N weighting matrices calculated for robot i are horizontally augmented to form W i (f b ). Using matrices to weight the relative unit vectors allows for more complicated behaviors than simple scalar weighting.
The weighted sum is premultiplied by K (θ b ) to determine the velocity for the behavior. Here, θ b is an angle, and K is the standard rotation matrix about the z-axis. A θ b = 0 deg will move the robot along the reference vector, while θ b = 180 deg will move the robot along the negative of the reference vector. Other angles of θ b are possible but are not used in this paper. Therefore, the b th behavior in an RPS composite behavior is specified by the weighting function and angle pair (f b , θ b ).
Reference [6] provides a detailed explanation of architecture and equations. It also establishes a library of common RPS behaviors found in literature, including several stable aggregation behaviors, a coverage behavior, a formation maintenance behavior, and several new orbiting behaviors. RPSs are reactive to their environment and the current swarm state. For this reason, the architecture design is very mathematical, avoiding complicated or nested conditionals to allow for easy scaling up and down in the number of robots. The architecture also uses simple functions and superpositions, avoiding long computational delays associated with consensus filtering [4]. The architecture uses only the current swarm state, unlike Particle Swarm Optimization (PSO). PSO saves the past local and global best solutions, which leads to non-reactive control, especially in rapidly changing environments [5].

III. ELLIPTICAL ATTRACTION ALGORITHM
RPS elliptical attraction is formed using two base behaviors: 1) an elliptical weighting function, f EQ , and 2) a nearest neighbor dispersion behavior, fñ. The following two subsections discuss each base behavior individually before describing the composite behavior.

A. WEIGHTING FUNCTION:f EQ
The characteristic equation of an ellipse is expressed in (4) for points (x, y) on an ellipse.
In (4), (x c , y c ) are the coordinates of the center of the ellipse, α is the angle between the semi-major axis and the x-axis, a is the semi-major axis, and b is the semi-minor axis. For any point (x, y), (5) calculates the parameter EQ.
An EQ of one indicates that the point is on the ellipse. If EQ is less than one, the point is inside the ellipse, whereas if EQ is greater than one, the point is outside the ellipse. The parameter max(EQ) indicates how well a swarm has aggregated to the commanded ellipse. max(EQ) is the maximum EQ value for a robot in the swarm. Since EQ indicates the position of a robot relative to the ellipse, max(EQ) provides a good indicator if the swarm has not fully spread through the commanded ellipse (max(EQ) < 1), has fully expanded to the ellipse (max(EQ) ≈ 1), or has not contracted to the desired ellipse (max(EQ) > 1). The eccentricity, e, of an ellipse is related to the semi-major and semi-minor axes as in (6).
Equations (7) and (8) describe the weighting function, f EQ using ellipse parameters a, b, and α.
M EQ = a sin 2 α b Here, d ij is the distance between robots i and j, N is the number of robots within the ACR of robot i. The center of the ellipse ⃗ x c = (x c , y c ) is the local center of mass (CM) based on robot i's view of the swarm, ⃗ x CM ,i . The Iverson bracket notation indicates that the weighting function is only active when the robot is outside the desired ellipse.
The matrix M EQ was derived using the gradient of (5) as shown in (9). The elliptical attraction function moves the robot parallel to the negative gradient of the elliptical field centered at ⃗ x c with a given eccentricity and orientation.
The weighting function fñ is a nearest n neighbor weighting function as presented in Reference [6]. Equation (10) calculates the reference vector as the weighted summation of the relative position vectors pointing from robot i to its n closest neighbors. The distance between robot i and its n th closest neighbor is d n . Using n = 1 has been shown to distribute robots through an area well [6].

C. ELLIPTICAL ATTRACT COMPOSITE BEHAVIOR
RPS elliptical attraction is formed using two base behaviors: 1) an elliptical attraction with (f EQ , 0 deg) with desired ellipse parameters, a, α, and b and 2) the nearest neighbor dispersion (fñ, 180 deg) with n = 1. The elliptical attraction behavior will bring the robots toward the desired ellipse. Once the robot is on or inside the ellipse, the elliptical attraction behavior turns off, and the nearest neighbor dispersion will distribute the robots throughout the ellipse. Intuition says that the scaling constant for elliptical attraction needs to be greater than that for the dispersion behavior to stabilize the composite behavior. Section VI-A discusses the selection of scaling constants in more detail.

D. LYAPUNOV STABILITY ANALYSIS
The stability proof of this elliptical attraction behavior is broken into three parts. First, the stability of a single generic behavior with the center of mass as an equilibrium point is proven. Second, the stability of a composite behavior incorporating the CM stable behavior is proven. Finally, these results are applied to the elliptical attraction behavior presented in this work to show that it is Lyapunov Stable.

1) CENTER OF MASS BEHAVIORS
Consider a single RPS base behavior of the form For simplicity, assume all robots are within d l of all other robots, so N = N and C i (d l ) = I 2N . So the behavior command velocity becomes VOLUME 11, 2023 Combining this with the assumed general form of f b , the command velocity simplifies to the following where ⃗ x CM is the center-of-mass position vector.
Equation (1) calculates the robot command velocity for a single RPS base behavior with no external inputs.
This is true for all robots in the swarm, indicating that ⃗ x CM is an equilibrium point. Next, let us look at the stability of this behavior. To begin, we use the coordinate transform in (16) to express the positive definite Lyapunov candidate function in (18).
Since⃗ x CM is independent of the summation, it can be moved outside the summation.
Substituting for⃗ x i , the time derivative becomeṡ Therefore, as long as M is a positive definite matrix, thenV is less than zero, implying that the equilibrium point, ⃗ x CM , is Lyapunov stable for this behavior.

2) COMPOSITE BEHAVIOR WITH CM STABLE BEHAVIOR
Next, we consider how this changes if the stable behavior is one of p base behaviors in an RPS composite behavior. There are also q external inputs. For simplicity, call the CM stable behavior b = 1, with a base behavior command velocity unit vector ofû 1i .
In the worst-case scenario, all other p + q − 1 behaviors are unstable, working against behavior b = 1. Then for the remaining p − 1 base behaviors with b ̸ = 1,û bi = −û 1i . Additionally, for all q external inputsû ui = −û 1i . The robot command velocity becomeṡ Applying the scaling constant constraint equation Now that we have derived an equation for⃗ x i , let us consider the equilibrium points. We know thatû 1i = ⃗ v 1i ||⃗ v 1i || . Combining this, (3), and (31), it is clear that since ⃗ x CM was an equilibrium point for the base behavior, it is also an equilibrium point of the composite behavior.
Using the same Lyapunov candidate function as before, Substituting in for the robot velocity, the time derivative of the Lyapunov candidate function becomeṡ We have already shown that M must be a positive definite matrix for behavior b = 1 to be Lyapunov stable. So forV to be less than zero, k 1 must be greater than 0.5.
In conclusion, a CM-stable RPS behavior causes a composite RPS behavior to be stable if and only if its corresponding scaling constant is greater than 0.5.

3) STABILITY OF f EQ
Now that we know that weighting functions of the form N are stable for positive definite matrices M , we must show that M EQ is positive definite. Since M EQ is a symmetric matrix and the eigenvalues b/a and a/b are both positive, M EQ is a positive definite matrix. Therefore, the proposed elliptical attraction composite behavior is Lyapunov stable if k b_EQ is greater than 0.5.

IV. ELLIPTICAL ATTRACTION ADAPTIVE SIZING
Elliptically attracting is useful, but a priori knowledge of the perfect size for the ellipse might be challenging. Instead, it would be helpful to adaptively size the ellipse based on the swarm state. This version of adaptive sizing scales the commanded semi-major axis while allowing the eccentricity and orientation of the ellipse to be specified as functions with respect to time.
Adaptive sizing can be done using the state machine in Fig. 2. The state machine uses a given minimum semi-major axis a min and an ACR, and the original commanded semi-major axis a 0 as well as functions f (x, t) and g(x, t). The function f (x, t) calculates the commanded ellipse orientation α while g(x, t) calculates the ellipse's eccentricity. Both f and g are functions of the swarm state and time.
As the state machine scales the semi-major axis, the value is clamped to be within a min and 3 8 d l . The lower bound a min is an engineering choice based on the application. The upper bound prevents the swarm from expanding out too far, losing communication between robots, and breaking the elliptical attraction behavior, as shown in Sec. VI-C1.
The state machine uses a function h that can be a function of variables such as swarm state or time. The ellipse should shrink if h is -1 and expand if h is 1. There are three transitions for the state machine based on h and the max(EQ) of the swarm. First, if h is -1, the swarm should only scale down if the swarm has reached the currently commanded ellipse indicated by having a (max(EQ) − 1) less than or equal to some tolerance, . When the robots have fully attracted to the ellipse, the obstacle avoidance or the nearest neighbor dispersion may push robots outside the ellipse for a time step.
should be selected to provide a buffer for when this happens. If these two conditions are satisfied, the new semi-major axis is a k+1 = a k − γ d a 0 . Second, if h is 1, then the swarm should scale up, with a k+1 = a k + γ u a 0 , only if it has max(EQ) greater than or equal to one. This condition prevents the swarm from scaling up on successive scaling up steps before the robots have fully expanded into the currently commanded ellipse. γ d and γ u are positive constants representing the fraction of the original semi-major axis the algorithm should scale by on each step. Third, if neither of the previous two cases is met, the state machine does not scale the semi-major axis.

V. CASE STUDY DESIGN
Next, a series of simulated and on-hardware case studies explore more details of the elliptical attract and adaptive sizing algorithms. All case studies were performed using Santa Clara University's Swarm Behavior Simulator (SBS) [22]. SBS is a MATLAB and Simulink-based platform capable of performing simulated and hardware-in-the-loop trials. The hardware trials connect to Santa Clara University's Decabot Testbed, a 3.5 m by 8.25 m workspace that localizes robots using an OptiTrack system [23]. The testbed commands up to 12 custom holonomic Decabots, as seen in Fig. 3. Fig. 4 shows the Decabots in a portion of the testbed. These holonomic robots allow for testing behaviors on hardware without having to deal with non-linear robot dynamics, as is customary with RPS. SBS also performs two types of simulations, single and batch simulations with randomized initial conditions, using simulated Decabot dynamics. These batch simulations allow for the comparison of swarm performance for different parameter sets, as seen later in this paper.
SBS uses an obstacle avoidance function independent of the RPS architecture to prevent collisions between robots. If two robots are within a set distance of each other, d OA , then (36) provides the robot command velocity. Otherwise, (1) calculates the robot command velocity.
In order to determine how well the robots attract to the commanded ellipse, two performance parameters are used 1) max(EQ), as discussed previously, and 2) the ellipse parameters and associated area for the final positions of the robots. The area and ellipse parameters are calculated using the Khachiyan Algorithm for prescribing a minimum volume enclosing ellipsoid [24]. The code to calculate the ellipsoid is adapted from the MATLAB File Exchange code in Ref. [25].   This paper uses two statistical tests to compare batches of simulations. The first method is the one-way ANOVA. The one-way ANOVA analyzes the difference between the means of more than two groups for one independent variable. The one-way ANOVA null hypothesis is that there is no difference among the group means. This test returns a p-value. If the p-value is less than the standard threshold of 0.05, the null hypothesis is rejected, indicating that at least one mean is different, and a post-hoc analysis is needed. The Tukey-Kramer post-hoc test will determine which groups have different means for both ANOVA methods. The Tukey-Kramer Test determines a p-value for each pair of groups. The means of the two groups are considered statistically different if p < 0.05. For all analyses, outliers were defined as 1.5 times the interquartile range above the 75 th percentile or below the 25 th percentile [26].

VI. GUIDELINES FOR SUCCESSFUL SWARM ELLIPTICAL AGGREGATION A. ADEQUATE SELECTION OF SCALING CONSTANTS
First, we look at selecting the scaling constants, k b_EQ and k b_ñ . The stability analysis in Section III-D showed that selecting k b_EQ greater than 0.5 should yield a stable behavior. Various numbers of robots were simulated attracting to an ellipse with an eccentricity of 0.9 and a semi-major axis of 10 m. The value of k b_EQ was varied from 0.2 to 0.9. The complete set of parameters is specified in Table 1. For a given number of robots and value of k b_EQ , 100 simulations were run, and a max(EQ) was determined for each simulation.
First, Fig. 5 shows the max(EQ) for each pair of number of robots and k b_EQ . The most prominent factor influencing max(EQ) is k b_EQ . To further investigate this, Fig. 6 shows the distribution of max(EQ) for each value of k b_EQ . Values of k b_EQ less than 0.5 result in large max(EQ) values due to the instability of the behavior. When k b_EQ is 0.5, max(EQ) is still relatively large compared to when k b_EQ is greater than 0.6. These results agree with the stability analysis in Section III-D.
Next, we examine for k b_ellipse > 0.5 whether increasing the scaling constant value yields better aggregation to the commanded ellipse. The mean max(EQ) for all distributions is less than 1.1, indicating good adherence to the ellipse. A one-way ANOVA analysis determined that no mean differed from the others (f (3) = 2.568, p = 0.0529). Therefore, increasing k b_EQ for pure elliptical attraction does not improve performance. Therefore, the remainder of the case studies use a k b_EQ value of 0.9.

B. DEMONSTRATION OF PERFORMANCE
The next series of simulations demonstrate that the commanded ellipse is reached for a range of angles, eccentricities, and semi-major axes. The parameters for the simulations are summarized in Table. 2. One hundred simulations were run for each set of angles, semi-major axes, and eccentricities. Simulations with commanded semi-major axes of 25 and 50 used runtimes of 500 and 700 seconds, respectively, to provide ample time for the robots to disperse through the ellipse. The robots started with random initial conditions in a 10 m by 10 m area for smaller semi-major axes and in a 40 m by 40 m area for semi-major axes of 25 and 50. In order to remove the division by zero when calculating the percent error when α = 0, the equivalent α of 180 deg is used.
First, Table. 3 calculates the mean and standard deviation for the distributions for each of the commanded semimajor axes. These results show good adherence between the commanded and calculated semi-major axes with a percent error (PE) of the mean less than 0.1% from the commanded.  The mean semi-major axis is expected to be lower than the commanded one because of the elliptical weighting function pushing robots in and the bounding ellipse calculating the tightest ellipse possible.
Second, Table 4 shows the mean and standard deviations for the five eccentricities. There is good adherence between the desired and calculated eccentricities with a percent error in the mean of at most 1%. The standard deviations and percent error of the mean eccentricity are more significant for smaller eccentricities because, for less eccentric shapes, there is more variation in the bounding ellipse used to calculate ellipse parameters due to the robot positions.
Third, let us examine how well the angle of the ellipse tracks the commanded ellipse. For less eccentric ellipses, we expect more error in the angle calculated by the bounding ellipse. Again, this is because the bounding ellipse calculates the tightest ellipse based on the individual robot positions but these low eccentricity ellipses are very circular, making calculating α more error-prone. For this reason, we determine the mean and standard deviations of the angles separated by commanded angle and eccentricity. Table 5 summarizes each distribution's mean and standard deviation as well as the percent error in the mean of each distribution. As expected, there were more significant standard deviations and larger percent errors for the calculated angles for less eccentric ellipses. However, overall the mean calculated angles matched the commanded angles to within a degree and with percent errors less than 1%, indicating good adherence to the commanded angle. Fig. 7 shows the distributions of max(EQ) for each semi-major axis and eccentricity. The color represents the mean value of max(EQ) with lighter shades indicating results that adhere less well to the commanded ellipse. The size of the circles represents the standard deviation in the values of max(EQ). Overall, there is good adherence to the ellipse with a maximum mean max(EQ) of 1.0067 and a maximum standard deviation of 0.0030. The smaller the semi-major axis and the more eccentric the shape, the higher the mean and SD of max(EQ). The reason for this is twofold. First, as the semi-major axis decreases and the eccentricity increases, the area of the commanded ellipse decreases, making it harder for all the robots to fit in the commanded ellipse while honoring VOLUME 11, 2023 TABLE 5. Mean, standard deviation, and percent error of the mean of calculated angle α for the commanded angle and eccentricity in degrees.

FIGURE 7.
For a given semi-major axis and eccentricity, the color of the represents the average max(EQ), and the size represents the standard deviation of max(EQ).
the obstacle avoidance range. Second, for these ellipses with smaller commanded areas, (5) gives larger values for the same disturbance outside the ellipse.
Overall, these results indicate that the proposed elliptical attraction function can attract swarms to a wide range of commanded ellipses. For simplicity, since we know the elliptical attraction function tracks the commanded parameters very well, we will often only use the final area and max(EQ) as performance parameters. The final area indicates how well the commanded semi-major axis and eccentricity were tracked, and max(EQ) indicates how well the robots cover the commanded area. Additionally, most subsequent analyses use α = 0 since this is just a coordinate rotation.

C. PREDICTING SWARM SUCCESS FROM COMMANDED PARAMETERS
Next, we investigate how to predict if a swarm can elliptically attract to a region, which depends upon the number of robots, the obstacle avoidance distance, the artificial communication range, and the area of the commanded ellipse.

1) INFLUENCE OF ACR
A swarm's artificial communication range influences its ability to attract to the desired ellipse successfully. For this elliptical attraction behavior, each robot calculates the local CM of the visible robots and uses that as the center of the commanded ellipse. If the ACR is too small, the robots cannot communicate with robots on the other side of the commanded ellipse, and the performance will be degraded, as seen in Fig. 8. Table 6 summarizes the parameters for the three simulation tests. Twenty robots are commanded from the same initial conditions to an ellipse with a semi-major axis of 10 m for three ACR ranges. An ACR of 5 m is too small, and the robots distribute through the space but not elliptically. When the ACR increases to 10 m, the resultant shape is more eccentric, but because the robots on either end can only see half of the ellipse, they do not converge to have the proper semi-major axis. Finally, when the ACR is 20 m, the robots can barely see the ellipse's other end, and the correct semi-major axis and eccentricity are achieved.
To remedy this issue, a good rule of thumb is to only command robots to ellipses with a semi-major axis of at most 3 8 d l . This means that if two robots were perfectly on opposite ends of the ellipse, they would be 3 4 d l away from each other, providing margin.

2) COMMANDED AREA VS. FINAL AREA
Given an obstacle avoidance range and the number of robots, it might be impossible to fit all the robots in the commanded area. Table 7 presents the simulation parameters for this case study. Fig. 9 shows how, for an obstacle avoidance range of 1 m, and different commanded areas, as the number of robots increases, eventually, the robots cannot all fit inside the commanded ellipse. Interestingly, when the robots fail to attract to the desired ellipse in Fig. 9, the final areas form a plane. This indicates that there could be a relationship between the minimum area a swarm can attract to, the number of robots, and the obstacle avoidance distance. Sec. VI-C3 investigates the relationship between these parameters.    commanded to an ellipse with a semi-major axis of 0.1 m and an eccentricity of 0.9 for multiple obstacle avoidance ranges. This ellipse area is too small for the robots to fit in. Instead, the robots will attempt to minimize the area they cover while honoring the obstacle avoidance range, providing  a reasonable estimate of the minimum area that a number of robots and d OA can attract. For each d OA and number of robots, 100 simulations were performed. A summary of the simulation parameters is found in Table. 8 Fig. 10 shows that for a given d OA , the slope relating the number of robots to the final area is relatively linear. Using a linear best fit, the slope (ρ) and offset (b) were determined for each d OA and are summarized in Table. 9. Fig. 11 shows that the density versus d OA is linear on a log-log scale. Eqn (37) express the linear best-fit equation for the robot density for a given obstacle avoidance range, ρ(d OA ).
The density is approximately proportional to the inverse square of the obstacle avoidance range, which agrees with intuition. Equation (37) does not have a bias because as the VOLUME 11, 2023 obstacle avoidance range goes to infinity, we expect the robot density to go to zero. Then the number of robots N est that can fit in some area A min can be estimated as in (38). The average bias from Table 9 was used to generalize the bias in (38).
In Fig. 10, for each d OA and number of robots, there is a range of areas that the swarms attract to. There is a larger spread for the final area for larger values of d OA and a higher number of robots. The minimum area that the estimate equations calculate is roughly in the middle of the area distributions due to the use of the linear-best fit. Therefore, the estimate equations will sometimes underestimate the area a swarm can attract to for a given d OA due to initial conditions and the constant velocity controller. Using the data from Fig. 10, the estimated area underestimated the final area by at most 10% for all cases. Therefore, a margin of 10% should be used when applying the estimate equations right at the boundary of possible ellipses. The estimate equation will also sometimes overestimate the minimum final area for some initial conditions. This is not a concern because if the swarm can fit in a smaller area for a given number and d OA , the swarm will easily achieve the larger area.

D. APPLICATION OF ESTIMATE EQUATION
There are three different ways the above equations can be used: 1) determine how many robots can fit given an area and d OA , 2) determine the largest obstacle avoidance that allows a specified number of robots to fit in a prescribed area, and 3) determine the minimum area a specific number of robots can fit in given d OA . The following three subsections perform   Table. 10 lists the simulations' full set of parameters. If we include the 10% margin for the area, we reduce the area to 180 m 2 so that any underestimation will still allow the swarm to fit in a 200 m 2 area. The robot density remains the same, so only 51.22 robots should fit in the commanded area. Fig. 12 shows the distributions for max(EQ) for each eccentricity and number of robots. The color of each circle represents the mean of max(EQ) while the radius represents the standard deviation of max(EQ). For less than 65 robots, the max(EQ) is close to 1, indicating that the desired ellipse was reached. Additionally, there is minimal variation in the standard deviation, indicating that eccentricity did not influence the swarm's ability to attract to the ellipse successfully. The solid and dashed horizontal lines indicate the number of robots possible from the estimate equation without and with margin, respectively. Both numbers of robots are within the region of successful attraction.
At 65 robots, the standard deviation increases, but the mean max(EQ) is still very close to 1, indicating that the swarm is right on the threshold of fitting in the ellipse. As the number of robots increases over 65, the standard deviations stay relatively constant. However, the mean max(EQ) steadily increases, indicating that the threshold has been passed and there are too many robots to fit in the ellipse. Overall, this case study demonstrates that (37) and (38) predict the robot density and number of robots well.

2) APPLICATION TWO: GIVEN A AND N, FIND d OA MAX
The next application of the estimate equations uses the area of the desired ellipse and the number of robots to determine the greatest obstacle avoidance range. Here 25 robots are  commanded to an area of 150 m 2 . The equations predict a density of 0.13 robots/m 2 and a maximum obstacle avoidance range of 2.76 m. Including the 10% margin on the area, the area is reduced to 135 m 2 . The estimated maximum density is then 0.15 robots/m 2 , and the maximum obstacle avoidance range is 2.62 m. The simulation parameters are listed in Table 10. Fig. 13 shows the resulting mean and standard deviations of max(EQ) for each parameter set. For d OA less than 3, the mean of max(EQ) is very close to 1 with little standard deviation. The max(EQ) means and standard deviations increased as the obstacle avoidance range increased over 3 m. The estimated maximum d OA without margin, the solid black line in Fig. 13, is right on the threshold of the robots fitting in the ellipse. Using the 10% margin provides some buffer for the estimate equation. Overall, the estimate equations predicted the maximum obstacle avoidance range very well.

3) APPLICATION THREE: GIVEN d OA AND N, FIND MINIMUM A
Next, we look at how well the equation predicts the minimum area the robots can fit in for a given obstacle avoidance range and a number of robots. For 15 robots with an obstacle avoidance range of 2 m, the equations predict a robot density of 0.2575 robots/m 2 and a minimum area of 39.34 square meters. This area is increased by 10% to 43.28 m 2 to accommodate the margin. Table 10 shows the complete set of simulation parameters. The simulation results are displayed in Fig. 14. The simulations have large standard deviations for the commanded areas at or below 42 square meters for max(EQ). Even though the mean max(EQ) for a commanded area of 42 m 2 is close to 1, the high standard deviation indicates that the swarm is struggling to attract to the correct area for all initial conditions. For commanded areas above 42 m 2 , there is a significant decrease in the standard deviation of the ellipse, especially at high eccentricities. Additionally, the mean max(EQ) becomes less than 1.1, indicating that the robots have converged to the desired ellipse.
The estimate equations predicted that the 15 robots with a d OA of 2 m could attract to an area of 39.34 m 2 and an area of 43.38 m 2 with margin. Based on the mean and standard deviation of max(EQ), the estimate equation underpredicted the area to which the swarm could aggregate. The 10% margin provided the extra room to ensure that the commanded area was sufficient to allow the swarm to aggregate successfully. In conclusion, the overall equations provided a good prediction of the minimum area.

VII. HARDWARE VERIFICATION OF ELLIPTICAL ATTRACTION BEHAVIOR
The following subsections use a series of hardware-in-theloop case studies to demonstrate the performance of the elliptical attraction behavior.

A. INDIVIDUAL PARAMETER COMMANDING
The previous case studies all used constant ellipse parameters. Additionally, the motion of the CM of the ellipse was not VOLUME 11, 2023  commanded, so the swarm aggregated to an arbitrary position. This case study demonstrates the attraction functions' ability to track the commanded ellipse parameters while undergoing translations, rotations, scaling, and shaping. For each transformation, the swarm will attract to the initial commanded ellipse parameters for 45 seconds. Then the attraction function will track 1) a CM trajectory for translation as in (39), 2) a trajectory in α for rotation as in (40), 3) a trajectory in a for scaling as in (41), and 4) changes in eccentricity for shaping as in (42). All other ellipse parameters will be held constant in each stage to demonstrate individual transformations. The translation behavior uses an external input command velocity in addition to the elliptical attract behavior. The trajectory for the CM is expressed in (39), and the external behavior command velocity is expressed in (43). The parameters for the individual ellipse parameter commanding hardware trial are summarized in Table 11.  Fig. 15a shows the max(EQ) for the trial. The swarm stays within the commanded ellipse after the initial transient for the robots to converge to the desired ellipse, even as the ellipse translates or the parameters change. Fig. 15b shows the disturbance magnitude from the commanded CM trajectory. The CM tracks the commanded trajectory well with less than a 1 cm disturbance.
Finally, Fig. 15c shows the percent error in the commanded ellipse parameters for each transformation. The swarm converges to the commanded ellipse after 20 seconds. Then as the commanded ellipse undergoes translation, rotations, and scaling in the subsequent three transformations, the swarm tracks the commanded parameters well, with at most a 3% error in all of the parameters. In the final phase, where the eccentricity changes, there is a spike in the percent error in all ellipse parameters. This is expected because when the robots are inside the commanded ellipse, their relative positions are controlled by the nearest neighbor dispersion algorithm, and it takes time for the robots to distribute through the commanded area. While distributing, the swarm does not track the commanded ellipse parameters. Additionally, max(EQ) for this translation does not spike, indicating that the swarm is still inside the commanded ellipse. Overall this shows good adherence to the desired ellipse parameters while the swarm undergoes individual transformations of translating, rotating, scaling, and shaping.

B. COMPOSITE PARAMETER COMMANDING
The previous subsection showed that this attraction function could handle individual translations, rotations, scalings, and shapings. This case study varies all three ellipse parameters simultaneously according to the Eqns. (44)-(46) while also translating as in (47).
The parameters for the trial are summarized in Table. 11 and Fig. 16 shows the results for the case study. The vertical line marks when the ellipse parameters start to change, and the translation begins. Fig. 16a shows that after the swarm settles to the commended ellipse, max(EQ) stays very close to 1, indicating good adherence to the commanded area. Fig. 16b shows that the swarm CM  tracks the commanded trajectory well, with about a 1 cm disturbance. Finally, Fig. 16c shows the percent error in the three ellipse parameters. The semi-major axis percent error after the swarm converges is less than 2%. The percent errors in the angle and eccentricity are less than 4% but grow with time. This is expected because, for lower eccentricities, the bounding ellipse that measures the ellipse parameters is less accurate due to the small number of robots used in the trial. Additionally, for more circular shapes, the angle calculation is more error-prone. The spike in the ellipse parameters percent error that was present during the individual shaping trial is not present here because the elliptical shaping weighting function, f EQ , is also influencing the swarm motion, not just the nearest-neighbor dispersion, fñ. Overall, a percent error of less than 4% for all ellipse parameters while also undergoing translations indicates that this elliptical attraction behavior performs very well.

C. APPLICATION OF ESTIMATE EQUATION
In Section VI-C3, a relationship between the minimum area of an ellipse, obstacle avoidance distance, and the number of robots was determined. Then, a series of simulations showed that these equations held. Now, the next three subsections demonstrate on hardware that these equations hold.

1) APPLICATION ONE: GIVEN A AND d OA , FIND MAX N
This case study predicts how many robots can fit in an area given an ACR. The trial parameters are listed in Table 12. The final ellipse parameters, the minimum distance between robots, and max(EQ) are listed in Table 13. The estimate equations calculate that 10.12 robots should fit into an ellipse with a d OA of 1.1 m and an area of 6.1216 m 2 (a = 1.5 m, e = 0.5). Including the 10% margin, only 9.60 robots should fit in the area.
Based on max(EQ) and the percent error in the commanded ellipse parameters, twelve robots did not fit, and eight robots did. The max(EQ) for ten robots was 1.24, which indicates that the robots did not all fit in the ellipse. While the original estimate equations predicted that ten robots should have fit in the area, the 10% margin dropped the estimate to below ten robots, so it is not surprising that ten robots did not fit. This shows, especially when using the 10% margin, that the estimate equations can accurately predict the number of robots that can fit in an area.

2) APPLICATION TWO: GIVEN A AND N, FIND d OA MAX
This case study tests the largest obstacle avoidance range, d OA , that allows N robots to fit in a specified area. The equations predict a maximum d OA of 1.5 m for 12 robots to fit within an ellipse area of 15.71 m 2 (a = 2.75 m, e = 0.75). The addition of 10% margin reduces d OA to at most 1.43 m. The trial parameters are listed in Table 12. Table 14 lists the final ellipse parameters, minimum distance between robots, and max(EQ) for each trial. According to max(EQ) and the final ellipse parameters, commanding an obstacle avoidance range of 1.4 m allows the robots to fit in the area with the correct ellipse parameters while commanding an obstacle avoidance range of 1.6 does   not. These results satisfy both the original estimate and the estimate with margin. This shows that the estimate equations predict the maximum obstacle avoidance range well.

3) APPLICATION THREE: GIVEN d OA AND N, FIND MINIMUM A
This case study determines the minimum area a certain number of robots can fit in given d OA . For the trial parameters in Table 12, (37) and (38) estimate that 12 robots can attract to a minimum area of 3.84 m 2 (a = 1.98 m, e = 0.95) given an avoidance range of 0.75 m. Using the 10% margin increases the minimum area to 4.23 m 2 (a = 2.08 m, e = 0.95) We command an area above and below the limit keeping the eccentricity constant. The final ellipse parameters, max(EQ), and minimum final distance between robots are listed in Table 15. As expected, commanding a semi-major axis of 2.25 m fits all robots, while commanding a semi-major axis of 1.75 m will not as based on max(EQ) and the ellipse parameters. These results satisfy both the original equation and the equation with the margin.

VIII. ADAPTIVE SIZING HARDWARE VERIFICATION
The previous hardware studies had the swarm follow trajectories for the ellipse parameters. Adaptive sizing allows  the size of the ellipse to be varied without needing to specify a trajectory. This section tests the proposed adaptive sizing algorithm on hardware using h as expressed in (48). For the first 45 seconds, if the swarm has aggregated to within a tolerance of the commanded ellipse for 5 seconds, the ellipse should scale down. After 45 seconds, the ellipse should scale up if the max(EQ) has been greater than 1 for 5 seconds. The parameters for the simulation are shown in Table. 16. With each scaling step, the swarm scales down by 0.15 and up by 0.2. The minimum semi-major axis, a min , was selected to 0.45 m because Eqns. (38) and (37) estimated a minimum  possible semi-major axis of 0.4 m. The ACR was selected to limit the semi-major axis to 1.25 m.
The commanded and measured semi-major axis and max(EQ) with respect to time are shown in Fig. 17. max(EQ) spikes every time that the adaptive sizing state machine changes the ellipse parameters, but then always settles close to one for at least the expected 5 seconds. The semi-major axis tracks the commanded value very well when scaling up and down. For the first 45 seconds, the commanded semi-major axis never goes below a min . When the semi-major axis is 0.55 m, max(EQ) settles close to but never goes below the 1.05 threshold. Therefore a min should have been selected with more margin.
For the last 45 seconds, the commanded semi-major axis never goes above 3 8 d l , which is 1.25 m. The percent error in semi-major axis, eccentricity, and angle α, are calculated only for times when the swarm has settled to the ellipse, 1 ≤ max(EQ) < 1.05. This is because the swarm can take on any shape as it expands or contracts to the desired ellipse. Table 17 presents the mean and max percent errors for the adaptive sizing trials. The semi-major axis tracks very well with an average error of less than 1%. The eccentricity and angle have average percent errors of less than 5%, which is very good. The bounding ellipse that calculates the ellipse parameters is much more sensitive to individual robot positions for small numbers of robots with low eccentricities and small semi-major axes. Additionally, SBS uses a constant velocity controller, which could cause robots to overshoot the desired ellipse for a timestep. Therefore it is not surprising that these spikes are present.

IX. CONCLUSION
The elliptical coverage problem for robotic swarms has several unique solutions. Successful swarm elliptical coverage should have complete coverage of an ellipse with commandable parameters, should be simple enough to both scale to a large number of robots and handle the addition and loss of robots naturally, and be able to adaptive sizing the commanded ellipse. However, no existing method has all three properties, and most methods are generally complicated. Additionally, no guidelines exist for selecting swarm and ellipse parameters to ensure successful elliptical aggregation. This work first presents a new elliptical attraction behavior that meets all three tenants of swarm elliptical coverage. The Lyapunov stable elliptical attraction behavior for Reactive Particle Swarms has good adherence to the commanded ellipse parameters and coverage of the area. Second, an adaptive sizing algorithm for the elliptical attraction behavior is presented. With little error, the dynamic control can rapidly adjust the swarm's position to track the changing parameters. Third, guidelines for selecting swarm parameters of the number of robots and communication limits, as well as all ellipse parameters, are presented to ensure the swarm can aggregate to the commanded ellipse. All results are validated with simulated and on-hardware case studies. Future work will use this elliptical attraction behavior and the associated swarm and ellipse selection guidelines to solve several navigation problems.