Fourier-Based Multi-Agent Formation Control to Track Evolving Closed Boundaries

The automatic monitoring/tracking of environmental boundaries by multi-agent systems is a fundamental problem that has many practical applications. In this paper, we address this problem with formation control techniques based on parame tric curves that represent the boundary’s feedback shape. For that, we approximate the curve with truncated Fourier series, whose finite coefficients are utilized to characterize the curve’s shape and to automatically distribute the agents along it. These feedback Fourier coefficients are exploited to design a new type of formation controller that drives the agents to form desired curves. A detailed stability analysis is provided for the proposed control methodology, considering both fixed and switching multi-agent topologies. The reported numerical simulation and experimental studies demonstrate the performance and feasibility of our new method to track closed boundaries of different shapes.

Formation controls can be categorized (based on the agents' sensing capabilities and interaction topologies) into the following three strategies [4]: position-based, displacement-based, and distance-based controls.Position-based strategies (which are the most widely used among the three) calculate the desired formation for the MAS using absolute position vectors [4].In this paper, we extend this basic idea and show how finite Fourier coefficients [5] (computed from a parametric contour defined in absolute coordinates) can be instead used to specify a desired formation, a strategy that effectively solves the lack of an explicit position target for each agent.
Formation controls are typically used to coordinate the motion of a variety of systems such as spacecraft [6], unmanned aerial vehicles (UAV) [7], [8], mobile robots [9], vessels and underwater systems [10], etc.One practically important application of formation control is the spatial monitoring of environmental boundaries, which is particularly critical e.g., to sample the spread of oil and dangerous chemicals during accidental spills [11], or to contain the boundaries of wildfires [12], see Fig. 1.As these regions may possibly be large in size, it is difficult to monitor them with a single autonomous robot; The use of MAS can significantly improve the efficiency of this task.Our goal in this paper is precisely to develop controls for these types of evolving curve-based formations.

A. Related Work
Formation control is one of the most popular control problems in MAS (see [4] for a comprehensive review).
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
State-of-the-art methods include [13], which proposes a position-based controller for vehicles subject to secondorder dynamics; [14] proposes a similar position-based strategy but for aerial vehicles.Representative examples of displacement-based controllers include [15] which presents a method to coordinate autonomous unmanned vehicles (AUVs) with optical sensors, and [16] which studies the maneuvering and robustness of undirected control topologies.Examples of distance-based strategies are presented in [17] and [18] for various conditions and applications.There are some works on formation tracking control, such as [9] and [19], yet, they mostly focus on driving followers to track leaders' trajectories rather than driving agents to track evolving parametric boundaries.
Formations computed with standard position-based strategies typically consist of discrete targets represented in absolute coordinates.Those computed with displacement/distancebased strategies do the control via relative inter-agent spatial configurations.However, these approaches may not be the most efficient way to distribute the agents along a boundary, whose most natural representation is a parametric curve.In our boundary-constrained formation problem, the agents must align along this parametric curve and automatically distribute according to the curve's arc length, an objective that clearly contrasts with that of traditional formation controls.Various related works have addressed similar contour-based scenarios.Some early examples include [20], which propose contour-based controls for homogeneous AUVs to track profiles of plumes.However, these methods require gradient measurements, which are difficult to obtain in practice.The method in [21] adopts the idea of active contour models from the image processing to compute a motion planner that guides robots into boundaries.The algorithm in [22] uses the advection-diffusion equation to model the dynamics of plumes and to control a team of robots that monitor them.The method in [23] describes a solution to estimate time-varying boundaries with multiple robots.More examples can be found in [24] and [25], where plume tracking and wildfire boundary monitoring with MAS are studied.Most previous works in the literature mainly focus on detecting/estimating the boundaries, do not exploit the properties of parametric curves in their design, have not considered the situation in which robots interact with switching topologies.Furthermore, they mostly focus on theoretical analysis, and very few conduct real-world experiments to verify the proposed methodology.The aim of this paper is to fill these gaps in the literature.

B. Our Contribution
In this paper, we propose a new control method to autonomously align a network of agents into an evolving parametric curve while keeping a desired arc length separation among them.The curve is represented with truncated Fourier series, whose "feedback coefficients" are used to establish a shape control loop.The main contributions of our work are summarized as follows: 1) A new representation of formation tasks with homogeneous MAS based on finite Fourier coefficients.
2) A new formation controller that explicitly uses Fourier coefficients to drive the agents towards the desired curve under fixed and switching topologies.3) A rigorous stability analysis, numerical simulation and experimental study to investigate the properties and validate the performance of the proposed method.Compared with existing works, our approach doesn't require explicit position targets for the agents.Instead, our proposed control method defines the formation in terms of a parametric curve (which can be simply obtained from an image) and drives the robots towards this target shape based on errors of Fourier coefficients (this represents a new type of servo-loop for formation control tasks).

C. Notation
Matrices and vectors are denoted as bold capital and small letters, respectively, e.g., M and m.We use [M] i j to denote the ith row jth column entry of a matrix M, and [m] i to denote the ith element of a vector.I denotes the identity matrix of appropriate dimension and 0 m×n is the m × n matrix of zeros.We use λ min (M) to denote the minimum eigenvalue of a matrix M.

D. Organization
The remainder of this paper is organized as follows.Section II models the parametric curves.Section III derives the proposed formation controller.Section IV presents the conducted simulation and experimental results.Section V provides final conclusions.

A. The Observed Boundary
Our interest in this paper is to study formation control strategies to drive agents toward closed planar curves.Therefore, the first step is to obtain the reference boundary.In actual missions, boundaries could be predefined or obtained from top-view images of a scene.Without loss of generality, we assume that these boundaries can be characterized as a time-dependent parametric equation of the form [26]: for s ∈ [0, 1] as the nondimensional parameter, and t as the time variable.c x (s, t) and c y (s, t) are the 2D position coordinates of a point along this curve assumed to be measurable.Since the curve is closed, it satisfies c(0, t) = c(1, t).

B. Curve Approximation
An intuitive way to automatically drive the agents and distribute them uniformly along the boundary of interest is to parameterize it by its arc length [27].However, a very accurate characterization of the observed contour (1) may require a large number of frequency terms, which invariably leads to complicated (and high dimensional) geometric representations.To deal with this issue, in our method, we use an approximation of the curve, which we denote as c for a new parameter where H is the number of harmonics, and a h , b h , c h , d h , e, f are the Fourier coefficients.The term c = [ cx , cy ] T denotes the 2D coordinates of a point along the approximated curve.The parameter s is defined as the normalized arc length of the observed curve (1), and it satisfies: for L(t) > 0 as the time-varying perimeter of the curve c, l(s, t) as the metric length of a curve section with the parameter range of [0, s].This approximation is depicted in Fig. 2.
For convenience in further processing, we rewrite (2) as a linear regression equation of the form: where ξ = [ρ 1 , . . ., ρ H , e, f ] T ∈ R (4H +2)×1 denotes the vector of parameters that characterize the curve, with ρ h = [a h , b h , c h , d h ] T (h = 1, . . ., H ) as the Fourier coefficients, and [e, f ] as the offset.The regression-like matrix G(s) = [g 1 , . . ., g H , I] ∈ R 2×(4H +2) contains the harmonic terms, for In (4), ξ is unknown and needs to be calculated.For that, we sample discrete points along the observed curve and stack them into a point sequence C ∈ R 2N ×1 of the form: where N > 0 is the total number of the sampled points.
For each point c p ( p = 1 . . .N ) in (6), we can obtain an arc length parameter sp by using (3), thus obtain a matrix G(s p ).We then stack these matrices into a similar sequence of the form: We calculate the Fourier coefficients by Note that the number of sample points along the curve must satisfy N > 2H + 1 to guarantee the existence of (G T h G h ) −1 , and thus, the boundedness of ξ .

III. FOURIER-BASED FORMATION CONTROL A. Graph Theory
The topology of the interactions among agents can be depicted by a graph, where agents are the nodes of the graph and interactions are the edges of the graph.The graph is denoted as G = (V, E, a i j ), where a i j are the weights of the interactions with the subscript i j indicating the i-th and j-th agents, V = {1, 2, . . ., n} is the set of node indexes, n is the amount of agents, and E = {(i, j) ∈ V × V : a i j ̸ = 0} denotes the set of edges.Given a graph, it is called an undirected graph if a i j = a ji , otherwise, it is a directed graph [4].
For a specific agent i in the multi-agent system, we denote its neighbors by N i = { j ∈ V : a i j ̸ = 0}.The degree of agent i is defined as the sum of weights of its neighbors, denoted by d i = j∈N i a i j .The interconnection of agents is established by the Laplacian matrix L ∈ R n×n , where a i j > 0 if j ∈ N i and a i j = 0 otherwise [28].

B. Controller Design
The agents considered in this work are assumed to conduct planar motions and have nonholonomic dynamics.We use [x T i , θ i ] T to describe the configuration of agent i, where x i = [x i , y i ] T and θ i denote the center position and the orientation of agent, respectively.The dynamic model of agent i is given by [29]: where T represents the control input, with v i and ω i as the linear and angular velocities of the agent, respectively.
A common practice for the control of nonholonomic dynamics is to use the input/output feedback linearization [30], where a virtual control point is introduced to change the coordinates to simplify the control objective and stability analysis.The virtual control point is given by [31]: for an arbitrary scalar distance ℓ > 0 that translates the agent's position to an arbitrary close location, see Fig. 3(a).
By using these new coordinates, we can construct the state vector xi = [ xi , ȳi ] T , whose time derivative yields a modified dynamic model of the form: Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.for a full-rank matrix defined as: By defining a new control input ūi = R i (θ i )u i , the position dynamics in ( 12) can be reduced to a single integrator ẋi = ūi .This enables us to represent the position dynamics of the multi-agent system as follows: for an extended position vector The control input of the original dynamics can be simply computed as To ensure that the agents have an equal arc length separation between neighboring agents along the curve, we uniformly discretize the arc length parameter s according to the number of agents, that is, the discretization step is ds = 1/n.Then, we obtain a sequence of arc length parameters corresponding to the distribution for each agent: where the ith length is calculated as si = (i − 1)ds.The assigning of arc length parameters is shown in Fig. 3(b).Substituting the sequence s into G(s), we similarly obtain the sequence of harmonic matrices corresponding to each agent as follows: Problem Statement: By using the extended position x and matrix (16), we can define the agents' position error: and the curve's coefficient error: where Ḡ+ denotes the pseudoinverse of Ḡ.The control objective is to design a Fourier-based feedback control input ū for the dynamics ( 14) that ensures the tracking of an evolving boundary c(t), i.e., the convergence of xe and ξ e as t → ∞, with arbitrary initial conditions.Remark 1: Note that the term Ḡξ can be interpreted as the target positions of the agents along a curve.In our proposed method, we codify the target shape in terms of frequency terms (i.e., the Fourier coefficients) ξ , which can be obtained e.g., from an image contour.
Remark 2: The rationale of ensuring the convergence of xi instead of x i .As is shown in Fig. 3, we can always ensure that the virtual control point is within the body of the agent by appropriately selecting the distance ℓ.In that case, the convergence of the virtual control point naturally means the agent's convergence to the tracked boundary.However, it should be noted that ℓ can not be zero, otherwise, the existence of the inverse matrix R −1 does not hold.
Since the interaction topology of the agents is represented by a graph G with the Laplacian L, by combining ( 14) with ( 16) we can design the following Fourier-based formation controller to drive the agents towards the reference curve: where k 1 , k 2 > 0 are feedback control gains.The first term on the right-hand side of ( 19) drives the agents to form the target shape of the curve.The remaining terms compensate for the relative distances between the agents and the reference curve.

C. Stability Analysis
We can rewrite the controller (19) and express it as the control input ū for the whole multi-agent system: where L = L ⊗ I is the extended Laplacian matrix, and ⊗ denotes the Kronecker product.The stability analysis of the proposed controller is based on the following assumptions: Assumption 1: The graph G that captures the interaction topology of the multi-agent system is strongly connected and undirected [2].
Assumption 3: As the observed curve evolves, the approximated curve is updated with a constant sampling period using the processes described in Sec.II-B.The observed curve evolves at a slow speed such that we can regard the coefficients ξ as constant during the sampling period.
With Assumption 1, we ensure that the Laplacian of the graph G is symmetric positive semi-definite, i.e. the eigenvalues of L are all non-negative.According to the property of the eigenvalues of a Kronecker product (see Theorem 4.2.12 in [32]), we can obtain the eigenvalues of L are all non-negative as well, i.e., L is positive semi-definite.
With Assumption 2, we know that when the number of agents n ≥ 2H +1, Ḡ is full column rank and its pseudoinverse is calculated by Ḡ+ = ( ḠT Ḡ) −1 ḠT .When the number of agents n < 2H + 1, Ḡ is full row rank and its pseudoinverse is calculated by Ḡ+ = ḠT ( Ḡ ḠT ) −1 [33].Based on Assumption 3, we separate the stability analysis into two parts.First, the stability of the controller is analyzed with time-invariant curves (i.e., set-point regulation).Then, the analysis is extended to tracking evolving curves.
According to the relationship between the number of agents and the number of harmonics, we separate the analysis with time-invariant curves into "many" agents case (n ≥ 2H + 1) and "few" agents case (n < 2H + 1).
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
Proposition 1 ("Many" Agents Case): Consider a MAS composed of n agents with dynamics ( 14) and a strongly connected interaction topology G. Given a time-invariant curve (4), the controller (20) ensures the asymptotic stability of the position and Fourier coefficients errors xe and ξ e , when the number of agents is n ≥ 2H +1 and the control gain k satisfies for = λ min L Ḡ Ḡ+ + Ḡ Ḡ+ L .Proof: The proof contains two steps.First, we prove the convergence of xe .By computing the derivative of xe we obtain the following relations: where , which is non-symmetric.Now, let us consider the following symmetric matrix: It can be seen that both items on the right side of ( 23) are symmetric matrices.According to the Theorem 4.3.1 in [34], we know that the eigenvalues of F 1 + F T 1 satisfy the bound: We can see that λ min (F 1 +F T 1 ) > 0 when k satisfies (21), which implies that all the eigenvalues of (F 1 + F T 1 )/2 are positive, thus, ensuring that F 1 is positive definite.The solution to the closed-loop dynamic system (22) with initial condition xe (0) can be computed as follows: which shows that the position error of the agents is asymptotically minimized, i.e., xe → 0 as t → ∞.
The second part analyzes the convergence of the Fourier coefficients error ξ e , which drives the curve-aligning actions in the controller (20).Note that since n ≥ 2H + 1, the pseudo-inverse of Ḡ is computed as Ḡ+ = ( ḠT Ḡ) −1 ḠT .By using ( 18) and ( 17), we obtain the following: As Ḡ+ has a full row rank, xe → 0 implies that ξ e → 0. Proposition 2 ("Few" Agents Case): Consider a MAS with n agents, dynamics (14) and strongly connected topology G. Given a time-invariant curve (4), the controller (20) ensures the asymptotic stability of xe and the boundedness of ξ e , when n < 2H + 1.
Proof: As Ḡ+ = ḠT ( Ḡ ḠT ) −1 , the closed-loop system satisfies the following: for a positive definite and symmetric matrix The solution of ( 27) is which proves the asymptotic stability of xe .To prove the boundedness of the Fourier coefficients error ξ e , we use (17) and ( 18) to obtain the relation: where by using (28), we can show that This expression implies the boundedness of ξ e , since ξ is constant for time-invariant curves.Remark 3: Proposition 1 proves that both xe and ξ e asymptotically converge to 0, which implies that ū → 0 as t → ∞ (i.e., zero motion at the equilibrium).Proposition 2 proves asymptotic stability of the linear position error xe but can only prove boundedness of ξ e .By substituting (30) into (20), we obtain Considering u = R −1 ū and since u i = [v i , ω i ] T , we can derive the following zero dynamics of orientation: which shows that the agents' orientation converges to bounded constant solutions.Remark 4: Note that static curves have constant Fourier coefficients.However, these coefficients are time-varying for evolving curves.Therefore, further analysis should be employed to analyze the stability of our controller in this situation.
Next, we show that the tracking position error vector xe is bounded on the basis of Propositions 1 and 2. We have shown the processes of approximation of the observed curve by the arc length at a specific time instant in Sec.II-B.The analysis of the tracking errors is based on the Assumption 3.
Before turning to the main analysis, we define several key terms useful for the theory's presentation.Let us denote the sampling period as t and the Fourier coefficients of the q-th sampling period as ξ q .We denote the difference of the Fourier coefficients of two neighboring sampling periods as ξ q , calculated by ξ q − ξ q+1 , express a sequence of coefficients differences with the following form: The maximum of ( 33) is denoted as ∥ Ḡ ξ ∥ max .To simplify notation, we introduce the matrix F ∈ {F 1 , F 2 }, and denote its minimum eigenvalue as ε = λ min ( F). Proposition 3 (Evolving Curves): Consider a MAS with n agents, dynamics (14) and strongly connected topology G. Given a time-varying curve, the controller (20) ensures the boundedness of the position tracking error ∥x e ∥ by the positive constant if the sampling period and the initial error satisfy the following condition: Proof: The detailed proof is given in Appendix A. Remark 5: Condition (35) points out that the sampling period can not be too small if we want to ensure the convergence of tracking errors.It can be derived from (34) that M can reach infinity if we select a very small t.Note that the existence of the natural logarithm in (35) requires ∥x e (0)∥ > ∥ Ḡ ξ ∥ max , which is generally easy to be guaranteed since we have made the slow motion assumption to the curve.

D. Switching Topology
In real-world applications, the communication among agents is usually limited by their relative distance from each other; This means that an agent can only communicate its state with agents within its local neighborhood.Network disruptions, such as unexpected hardware errors and poor signal intensity, can also cause communication failures.In this situation, the multi-agent system is said to interact through a switching topology.The aim of this section is to show that the proposed controller (20) is capable of handling switching topologies.
The presented analysis depends on the extension of Assumptions 1 and 2 in Sec.III-C.When communication failures occur, agents are divided into several subgroups where local communication is preserved.For that, our method relies on the following assumption: Assumption 4: All the interaction topologies within the subgroups that have more than one agent can be represented by strongly connected and undirected graphs.
We denote the set of all possible graphs as Note that the set G is finite because the number of agents is also finite.Proposition 4: Considering Assumption 4, the conclusions from Propositions 1 and 2 still hold for a multi-agent system interacting through switching topologies.
Proof: With a switch topology, the matrix L is not necessarily the Kronecker product of the Laplacian and identity matrices.However, we can always transform L into the following form: where L 1 , . . ., L α denote the Laplacians of the interaction graphs of the subgroups, and α is the number of subgroups that contain more than one agent.The identity matrix in Lt (i.e., the last term in the block diagonal) represents the isolated agents that do not interact with other agents.
From the aforementioned assumptions, we know that the matrices L 1 , . . ., L α are symmetric and positive semi-definite, therefore, Lt is symmetric and positive semi-definite.This condition on Lt enables us to directly follow the proofs of Propositions 1 and 2.
Proposition 5: Considering Assumption 3, a multi-agent system with switching topology and driven by the controller (20) can track an evolving curve with bounded position tracking errors.
Proof: According to the proof of Proposition 4, we can find a symmetric and positive semi-definite Lt for a switching topology.Accordingly, we can also obtain a positive definite matrix similar to F 1 or F 2 .Denote the matrix as F 3 , then we can define the minimum eigenvalues ε = λ min (F 3 ).Follow the procedure of the proof of Proposition 3, then we can prove Proposition 5.

A. Simulation Setup
1) The Observed Curve: We select a curve that evolves in terms of both its shape and the location of its center.We represent the observed curve c(s, t) by the following time-varying parametric equations: The curve starts as a circle and gradually evolves into a parallelogram-like shape.
2) Calculation of the Curve: We sample the points along the curve on a constant parameter step basis.Specifically, we set the constant parameter step as δs = 0.001, therefore, we collect N = 1001 sampled points The metric length of each curve section can be computed by the following curve integral [35] l(s p , t) = where s p = ( p − 1)δs.Note that the analytical computation of (39) usually is difficult.Therefore, in this study, we numerically compute the curve lengths by approximating (39) to: To approximate the observed curve, we set the number of harmonics to H = 6, which provides a compact shape feature vector ξ .
3) Formation Control: We perform two types of control tests.In the first, the agents interact with each other through a fixed topology.The interaction topology is established following the rules that the agent i connects to the agent i − 2, i − 1, i + 1, and i + 2. The weights of all interactions are set Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.to a i j = 1.Therefore, the Laplacian of the MAS is: otherwise. (41) In the second type of test, the agents interact with each other through a switching topology.The interaction topology switches among three modes: (1) all agents interact with their neighbors; (2) some agents interact with each other; (3) none of the agents interact with each other.
The duration of the control test is set to 200 seconds, the sampling period is t = 0.1 seconds, and the feedback gain of the controller is empirically set to k 1 = 2 and k 2 = 1.We set the distance ℓ as 0.01 m.The agents' positions and orientations at the beginning of each test are set as random values.
In the following sections, we present and analyze some representative simulation results to illustrate the validity of the proposed method.The readers are also referred to the supplementary video (https://vimeo.com/809466573)to check more simulation results.

B. Fixed Topology 1) n ≥ 2H+1:
We set the number of agents to n = 15 to ensure n ≥ 2H +1 in this simulation.The interaction topology of the multi-agent system is shown in Fig. 4(a).It can be calculated that ≈ −2.6822 × 10 −15 , therefore, the selected gains satisfy the condition in (21).The minimum sampling period is 3.4364 × 10 −4 sec.computed by (35).Therefore, the selected sampling period for the simulation is feasible.
The formation results and the evolution of errors are shown in Fig. 5.As shown in Fig. 5, the approximated curve obtained by (8) fits the observed curve well during the simulation; The agents successfully form the desired curve shape and track the evolving curve under the effects of the proposed formation controller.We record the norm of coefficient errors ∥ξ e ∥ in Fig. 5(d) and the norm of position errors ∥x e ∥ in Fig. 5

(e).
As is shown, the errors rapidly converge at the beginning stage and remain bounded at the stable stage.We also record the evolution of agents' orientations in Fig. 5(f).It can be seen that the orientations of all agents are bounded when the simulation reaches the stable stage, which meets the conclusion of Remark 3. The results in both figures follow Propositions 1 and 3, which illustrate the validity of the proposed controller.2) n < 2H+1: We set the number of agents to n = 8 to ensure n < 2H + 1 in this simulation.The interaction topology of the multi-agent system is shown in Fig. 4(b).It can be calculated by (35) that the minimum sampling period is 3.4013 × 10 −4 sec.Therefore, the selected sampling period for the simulation is feasible.
The formation control results are shown in Fig. 6(a)-(c).We can see that the agents successfully form the desired shape and track the evolving curve under the effects of the proposed controller.The evolution of errors and orientations are shown in Fig. 6(d)-(f).As shown in Proposition 2, the coefficient errors ξ e converge to a limit with the form of (30).For that, we define the following coefficient error metric whose evolution during the task can be seen in Fig. 6(d).The metric ξ m converges to zero as the simulation runs, which agrees with the conclusion of Proposition 2. The norm of position errors shown in Fig. 6(e) converge at the beginning stage and remain bounded as the simulation runs, as proved in Proposition 3. The orientations of agents are shown in Fig. 6(f).We can see that all agents have bounded orientations at the stable stage, which follows the conclusion of Remark 3.

C. Switching Topology
We set the number of agents to n = 15.The interaction topologies of the multi-agent system are shown in Fig. 7.As mentioned in Sec.IV-A.3, the interaction topology switches among three modes during the simulation.The agents interact

D. Experimental Validation
1) Setup: In addition to numerical simulations, we conducted some experiments to verify the performance of the proposed method in a real-world task.These experiments are conducted with the platform shown in Fig. 9, which is composed of nine (n = 9) Mona robots [36] with nonholonomic dynamics, a top-view camera to monitor the robots and obtain the feedback pose of the robots, a control PC to collect feedback data and send the motion commands.Wi-Fi communication (at a rate of 10 Hz) between the control PC and robots is built based on ROS.We conduct all experiments in an arena with a size of 1.8m × 0.8m.
In all experiments, the robots interact through a fixed topology, which is established following the same rule presented in Sec.IV-A.3.We empirically set the controller's gains to k 1 = 2 3 and k 2 = 2. Various representative results are presented in the following sections, yet, more experiments can be found in the supplementary video at https://vimeo.com/809466573.
2) n ≥ 2H + 1: In this experiment, we drive the robots to track an evolving ellipse whose size is growing over time.The equation of the observed curve is: At the beginning of the experiments, all robots have random initial positions and orientations.We set the number of harmonics in the curve approximation to H = 3 so that n ≥ 2H + 1 is ensured.It can be calculated that = −8.8818e×10−16 .Therefore, the selected control gains satisfy the condition in (21).As the minimum sampling period is 1.8 × 10 −3 seconds (which is computed by ( 35)), the selected communication frequency is appropriate for the system.Fig. 10 shows the results of this experiment.Fig. 10(a) depicts various screenshots of the robots during their motion, where we can see how they start from random initial conditions and gradually form the desired evolving shape while keeping track of it.The evolution of the coefficients errors, position errors, and orientations are shown in Fig. 10(b)-(d).From these plots, we can see the asymptotic properties and boundedness of the errors, which follow Proposition 1 and 3. Orientations also reach some bounded values at the stable stage of the experiment, which follows the conclusion of Remark 3.
3) n < 2H + 1: In this experiment, we drive the robots to track an evolving curve similar to (38).The curve starts from a circle, then gradually evolves into a parallelogram-like shape.The equation of the observed curve is: c x (s, t) = 0.0015[(0.5t* sin 4πs +200) cos 2πs +3t +300] c y (s, t) = 0.0015[(0.5t* cos 4πs +200) sin 2πs +300] (44) In this case, the number of harmonics is set to H = 6 to ensure n < 2H + 1.At the beginning of the experiment, all robots have random initial positions and orientations.As the minimum sampling period is 8.2674 × 10 −4 seconds (which is computed by ( 35)), the selected communication frequency is appropriate for the system.Fig. 11(a) depicts various screenshots of the robots during the experiment.The evolution of the coefficients error, position errors, and orientations are in Fig. 11(b)-(d).These plots similarly show the asymptotic properties and boundedness of the errors.We can also see that all agents have bounded orientations at the stable stage of the experiment, which agrees with the conclusion from Remark 3.

E. Discussion
1) Slow Motion: In the previous analysis, we have made the slow-moving assumption to the curve.In this section, we discuss the influence brought by this assumption.For that, we introduce a speed ratio to evaluate the moving speed of the curve, defined as r = v c /v max , where v c denotes the moving speed of the curve and v max denotes the maximum Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.linear velocity of the agent.It can be derived that r = 0 if we don't limit the linear velocity of agents.
Quantitatively measuring the moving speed of a curve is usually difficult.Therefore, instead of adjusting the curve, we set different limitations on the linear velocity of agents to examine the influence of different speed ratios.The results are shown in Table I.The settings of the simulations are the same as Sec.IV-B.1 except for the agents' maximum linear velocity and that agents set off from the same straight line with the same initial orientations.We set five different v max for the agents and record the norm of the position errors ∥x e ∥ at five different time instances.Table I shows the convergence speed of ∥x e ∥ decreases as v max also decreases.When v max = 0.1 m/s, the duration of the simulation is not enough for the MAS to reach the curve.Fig. 12 shows the slow convergence of ∥x e ∥ with v max = 0.1 m/s.

V. CONCLUSION
In this paper, we study the problem of tracking evolving parametric curves with multiple homogeneous agents.We approximate the observed curve as a Fourier series with arc length parameters, which ensures uniform inter-agent arc length separations.Then, a Fourier-based formation controller is designed to drive the agents to form and track the desired curves.We present a detailed stability analysis of the proposed formation controller and conduct simulations under fixed and switching topologies.Simulation and experiment results show that the proposed controller can drive the agents to form and track the desired curve under both fixed and switching topologies.We also discuss the influence of initial conditions and agents' maximum linear velocity on the performance of the proposed method.
However, the proposed method has various limitations.On the one hand, we have not considered collision avoidance between agents in this paper, which is important in real-world practice.On the other hand, we have only studied agents with nonholonomic dynamics.The formation control of agents with more complex dynamics, such as general non-linear dynamics or dynamics with disturbance [37], needs to be further studied.Future work also includes the formation control with different topology constraints, such as joint connected graphs and directed graphs.Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
Note that x1 e ( t) is obtained based on the truth that we regard ξ 1 as constant during the sampling period.Therefore, there is a difference between x1 e ( t) and the actual error.Since the end of the first period is the beginning of the second period, we know that the actual Fourier coefficients at the end of the first period are ξ 2 .Note that xe = x − Ḡξ , then we obtain the actual errors at the end of the first period: Accordingly, we can derive the actual error at the end of the q-th period: We can see that, as q → ∞ (i.e., t → ∞),

Fig. 1 .
Fig. 1.Conceptual representation of a MAS approaching a closed boundary.

Fig. 2 .
Fig. 2. The conceptual process of the curve approximation.
s ∈ [0, 1].The approximation c is obtained by fitting c with the following truncated Fourier series. cx

Fig. 9 .
Fig. 9.The setup of experiments: (a) the configuration of the platform; (b) the Mona robot.
e −ε t ∥ Ḡ ξ ∥ max .−ε t ∥ Ḡ ξ ∥ max , then we can derive that there always exists a positive constant M such that∥ a xq e ∥ = ∥x − Ḡξ q+1 ∥ ≤ M(52)holds as t → ∞, that is, the position errors xe are always bounded as t → ∞.Considering the condition(35), we can derive that M ≤ ∥x e (0)∥ always holds, which ensures the decrease of the errors.

TABLE I EVOLUTION
OF ∥x e ∥ UNDER DIFFERENT SPEED RATIOS