{'fraction_non_alphanumeric': 0.038929897973576565, 'fraction_numerical': 0.011001927688184683, 'mean_word_length': 4.445468509984639, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Optimal Time-Invariant Formation Tracking for a Second-Order Multi-Agent System

Given a multi-agent linear system, we formalize and solve a trajectory optimization problem that encapsulates trajectory tracking, distance-based formation control and input energy minimization. To this end, a numerical projection operator Newton's method is developed to find a solution by the minimization of a cost functional able to capture all these different tasks. To stabilize the formation, a particular potential function has been designed, allowing to obtain specified geometrical configurations while the barycenter position and velocity of the system follows a desired trajectory.


I. INTRODUCTION
During the last decade, designing and controlling groups of robots to achieve specific collective goals has drawn a considerable interest. An individual agent may be programmed to be fully autonomous, but because of physical and resource constraints, its abilities may be limited. On the other hand, groups of individuals exchanging information and optimally self-organizing may have larger capabilities. In nature, examples of interacting "swarms" abound [1]- [2] and, such inspired, multi-agent systems have been widely designed to be used in applications like vehicle coordination, exploration and mapping of unknown environments, cooperative transportation, surveillance and monitoring of dynamic scenes and crowds.
Related work: In the scientific community, a significant research effort has focused on the control of multi-agent systems and its mathematical foundation on graph-based networks [3]. Open challenges can be categorized as either formation control problems (with applications to mobile robots, autonomous underwater vehicles, satellites and spacecraft systems, automated highway and unmanned air vehicles, e.g. [4]- [5]), or other cooperative control problems such as role assignment, payload transport, air traffic control, timing, search and synchronization (e.g. [6]). Formation Control generally aims to drive multiple agents to achieve prescribed constraints on their states. Depending on the sensing and the interaction capabilities of the agents, a variety of formation control problems can be found in the literature (see [7] for an excellent overview on the subject). Nonetheless, creation and maintenance of a formation needs to trade off with single dynamics behaviors, collision avoidance, and dispersion [8]. In [9], an agent model that implements simple laws to abide these targets is introduced, thus determining the behavior of each agent within the group. These and similar laws often resort to a potential that determines the pairwise interactions among the agents, and along this line, many formation control algorithms adopt a formulation that contains an attractive part to maintain the swarm cohesion and a repulsive one to avoid agent collision (as the Morse potential in [10]- [11]). More generally, authors in [12] discuss results on Formation Control discerning among position-based, displacement-based and distance-based according to the types of sensed and controlled variables. In particular, the last approach is that adopted in this work and can be defined in a framework where inter-agent distances are actively controlled to achieve the desired formation, which is defined through the specific inter-agent distances. Finally, in this respect, it is worth to highlight that several novel achievements have been developing by using distance-based control laws and an extensive part of that research has been dedicated to Formation Tracking [13]- [14].
Contributions and outline of the paper: Many recent works developed for autonomous vehicles exploit Trajectory Optimization to perform maneuver regulation and motion planning, even in constrained environments. Practically, the employment of direct methods for the minimization of a cost functional represent a fundamental approach to provide solutions to this kind of problems: numerical tools, such as the PRojection Operator based Newton's method for Trajectory Optimization (PRONTO), have been successfully implemented in several instances [15]- [16]. At the light of this consideration, we propose here a novel application of PRONTO that aims at computing a solution for an optimization problem combining aspects of trajectory tracking, formation control and energy minimization.
In the remainder of the paper, Sections II and III describe how the specific problem of Formation Tracking can be formulated and solved by using PRONTO. Then, Sections IV and V show some numerical simulations, providing a validation of the proposed approach, interesting remarks and future directions for this research.

II. PROBLEM SETUP
In this section, we discuss the assumptions regarding the agents' dynamics and how these lead to a simplified analysis, formalization and implementation of a centralized algorithm to solve an optimal control problem that encapsulates distinct tasks such as trajectory tracking, time-invariant formation control and input energy minimization. In particular, we address this as the Optimal Time-Invariant Formation Tracking (OIFT) problem, with the aim of finding a potentialbased solution by minimizing a global cost functional able to capture all the different assignments simultaneously. Besides, it is crucial to highlight that communication topology constraints are not taken into account, since this paper lays the groundwork for a new optimization approach based on PRONTO in the field of Formation Tracking.

A. Agents dynamics
We assume that n > 1 robotic agents have been already deployed on a M -dimensional space, with M ∈ {1, 2, 3}. We also suppose that each agent i, for i = 1, ..., n, is aware of its absolute position p i ∈ R M and velocityṗ i ∈ R M and can be driven by regulating its absolute accelerationp i ∈ R M . For the sake of simplicity, a linear dynamics for the robots is adopted. With N := nM , the expressions of the state x ∈ R 2N and the input u ∈ R N of this linear system of mobile elements are given respectively by We finally assume that the state information is globally available to all agents, so that an estimate of the barycenter p i and velocityṗ B be available to each agent at each time instant. We set x B = p Bṗ B , with x B ∈ R 2M , to be an output for this linear system. Since we wish to govern positions and velocities of the agents by controlling their accelerations, a double integrator model is a suitable choice for this purpose: the second-order dynamics that follows can be represented via the linear state space State matrix A ∈ R 2N ×2N , input matrix B ∈ R 2N ×N and output matrix C ∈ R 2M ×2N in (1) are given by where I denotes the -dimensional identity matrix, Z 1× 2 indicates null matrices of dimension 1 × 2 and if 1 = 2 = the convention Z × = Z is adopted instead.

B. Cost functional minimization
Tracking a desired path x B,des = p B,desṗ B,des with the centroid x B while minimizing the energy spent by the input u and attaining a specified formation during this optimal control procedure is what we aim at. Let T be the trajectory manifold of (1), such that ξ = (x(·), u(·)) ∈ T . The general optimization problem can be formulated as follows: find a solution ξ such that represents the cost functional to be minimized in order to fulfill the three objectives previously stated. In (2), two different terms explicitly appear: the instantaneous cost and the final cost in which, matrix P 1 is set to be null in order to focus our analysis on the instantaneus cost (3) and simplify the overall framework. Each term in (3) is minimized with the aim to obtain the achievement of a specific task. Indeed, setting r ij = p i − p j 2 as the inter-agent Euclidean distance, we can define each contribution as follows: for the tracking task, for the input energy task and, given a family of potential for the formation task. Symmetric positive semidefinite matrix Q B ∈ R 2M ×2M , symmetric positive definite matrix R ∈ R N ×N and k F > 0 are constant weights that can be tuned according to given specifications in order to penalize the trajectory tracking, the energy spent by the inputs and the convergence to a desired shape respectively. Furthermore, with regard to potential functions σ dij in formation cost (6), each d ij represents the desired inter-agent distance between the pair of agents (i, j): an accurate choice of σ dij allows the usage of the optimization framework PRONTO, as explained next.

C. Potential based formations
As for the formation objective, we would like the system of agents to achieve desired linear, 2D or 3D shapes induced by a set of constraints of the form therefore, potential functions to accomplish this task are introduced. Let s ij = r 2 ij be the squared inter-agent distance between agents i and j. Among various possible choices, the potential function has been selected since it is polynomial in r ij with relatively low degree and differentiable with respect to s ij ≥ 0 until the second order. In particular, (8) is one of the C 2 functions with lowest polynomial order that can be adjusted by tuning constants k r > 0 and k a > 0 independently. These two hyper-parameters are purposely designed to be directly proportional to the repulsion and attraction action between agents respectively, playing a key role in the intensity regulation of the potential itself (see Fig. 1 in Sec. III-C).
The use of functions with the same properties of (8) can be justified by the fact that they can lead formations to satisfy the maximum number of feasible 1 relations in (7), since the dynamics is purposely steered to minimumpotential trajectories. Indeed, it can be easily proven that σ dij (s ij ) is nonnegative for all s ij and exhibits a unique global minimum point s ij , σ dij (s ij ) = d 2 ij , 0 , meaning that this potential vanishes when the desired distance d ij between the pair (i, j) is attained.

III. NUMERICAL METHODOLOGIES FOR OIFT
We provide here a brief introduction to PRONTO 2 , the numerical tool used for trajectory optimization, and clarify its adaptation to solve the OIFT problem introduced before.

A. Basics of PRONTO
In its simplest form, PRONTO is an iterative numerical algorithm for the minimization, via a Newton method, of a general cost functional (2). To work on a trajectory manifold T , one projects state-control curves in the ambient Banach space onto T by using a linear time-varying trajectory tracking controller. To this end, a nonlinear feedback systeṁ where matrix K acts as a time-varying controller, defines a nonlinear operator (C r when f ∈ C r ) P (ξ) : ξ = (α(·), µ(·)) → η = (x(·), u(·)) mapping bounded curves (in its domain) to trajectories.
Among the many properties of P(·), the most relevant in terms of the practical use of PRONTO refers to the fact that it makes possible to easily switch between constrained and unconstrained minimization to solve the initial problem, as where g(ξ) = h (P(ξ)).
Let us denote the Fréchet derivative of P(·) with the continuous linear operator such that the approximation P(ξ + ζ) ≈ P(ξ) + DP(ξ) · ζ holds for all ξ ∈ T . Then, PRONTO algorithm can be illustrated as follows: Algorithm 1 Projection operator Newton's method [19] 1: given initial trajectory ξ 0 ∈ T 2: for k = 0, 1, 2, ... do 3: Redesign feedback K(·) for P(·), if desired/needed 4: ζ k = arg min Step 4 in Alg. 1 encodes the search direction problem and can be solved by finding a solution to the time-varying linear quadratic (LQ) optimal control problem of the form where vectors a, b and matrices Q o , S o , R o ,Ā,B are time-varying parameters descending from the instantaneous cost l and the dynamics adopted. Moreover, vector r 1 and matrix P 1 are given constants depending on the final cost (4). When problem (9) possesses a unique minimizing trajectory, it can be solved resorting to the following differential Riccati equation in the matrix variable P: B. Application of PRONTO to the OIFT problem As previously mentioned, the assumptions made in Sec. II allow to obtain significant simplifications in the expression of many numerical quantities while preserving the possibility to perform interesting trajectory explorations for the system of agents. Now, we clarify how PRONTO variables and parameters declared in Alg. 1 have to be practically assigned and computed when this particular framework is adopted. Indeed, setting f (x, u) = Ax + Bu, and denoting the zero vector of dimension with 0 , the following list of equalities is immediately available: where symbols ∇ * and H * * indicate standard gradient and Hessian operators. Hence, a constant proportional derivative (PD) controller as is the simplest and most efficient choice to opt for, meaning that step 3 in Alg. 1 is not needed to be processed. In addition, let us examine the gradient and Hessian matrix of l f o d w.r.t. the state x appearing in (11) and (12) respectively. Their formal expressions are yielded by and, for all i = j, Note that the choice of having σ dij ∈ C 2 is required to use PRONTO. Furthermore, a nice expression for the diagonal blocks in the Hessian H pp l f o d can be provided by

C. Robust heuristic for PRONTO
As a matter of fact, the LQ problem in (10) requires a positive semidefinite matrix Q o to be solved; therefore, care has to be taken while search directions in Alg. 1 are computed. In reality, it can be noticed that expression (12) does not always satisfy this condition because of the general undetermined definiteness of the Hessian term. To this purpose, we decide to implement a suitable safe version of Q o , namely Q saf e o , by using a Gershgorin-circle-theorem-based heuristic that can be justified by the fact that potential functions with the same repulsive-attractive behavior of (8) never guarantee the positive semidefiniteness of matrix H pp l f o d . Indeed, the latter term is badly affected by the evidence that matrix Π ij is rank-1 for all the pairs (i, j) and the possibility of σ dij (s ij ) to be negative. As illustrated in Fig. 1, the qualitative graphic of (8) and its derivatives highlights that where the equalities in 1)-3) hold if and only if s ij = d 2 ij . The previous observations allow to state that, whenever two agents (i, j) are repelling each other, the Hessian H pp l f o d is not positive semidefinite. This implication leads us to neglect the variation of potential (8) w.r.t. the current distance for 0 ≤ s ij < d 2 ij to aim at searching a well defined descent is not required or σ dij (s ij ) > 0 then 6: is imposed, which, in most of the actual scenarios, becomes an infeasible 3 requirement to be completely verified for each pair (i, j). Set of constraints in (14) allows to explore specific formations, called d-configurations 4 , where potentials can be interestingly tested. Thus, if not differently specified, we adopt (14) as an assumption with d = 5 m during the dissertation. Whenever (14) becomes unfeasible, we relax the constraint to We have also chosen to keep the input weighting matrix R = r a I N , with r a = 1 m −2 s 4 , and the integration time T = 20 s fixed, as well as the structure of the output weighting matrix where q p , q v are nonnegative weights. Moreover, we have set the following parameters: k r = 10 2 , k a = 1, in (8) and k F = 10 −1 , in (6). In each numerical simulation presented, we have established to stop the execution of Alg. 1 by taking two actions: an interruption occurs when a preset maximum number of iterations M axIter > 0 is exceeded; the algorithm can terminate at the iteration k ≤ M axIter if, given a threshold ε > 0, the inequality −Dg (ξ k ) · ζ k < ε is satisfied. In particular, we have set ε = 10 −8 . Finally, constant gains for the PD controller in (13) have been assigned exploiting the linear dynamics of (1), as

A. Validity test for the developed algorithm
In the two simulations depicted in Fig. 2, the tracking of a straight line at the constant velocity v = 1 m s −1 is required to be performed by two formations composed of n = 3 and n = 4 agents in a planar (M = 2) and three-dimensional (M = 3) scenarios respectively. Denoting with vec(·) the vectorization operator 5 , the initial conditions are set as for the scenario in Fig. 2(a); for the scenario in Fig. 2(c). In addition, q p = 10 m −2 , q v = 1 m −2 s 2 and M axIter = 50 are assigned. It is worth to note that the functional g maintains strictly decreasing while the execution is running and both simulations correctly terminate before the maximum number of iterations is reached. Moreover, it can be numerically shown that instantaneous costs l tr in (5) and l f o d in (6) tend to zero as the time t goes to T , since d-configurations that satisfy all the constraints in (14) are attained and their centroids converge towards the desired trajectory.

B. Equilibria of the final configurations
In the group of simulations illustrated in Fig. 3, the equilibria of few particular formations are depicted: in this setup, agents are randomly deployed near the origin at the beginning and then asymptotically driven towards it. Therefore, q p = 10 m −2 , q v = 1 m −2 s 2 and M axIter = 100 are assigned. Assuming that a constraint of set (14) is satisfied if (15) holds true, it is important to note that, for some 5 By definition, this operator stacks the vectors assigned to its argument, i.e. vec w 1 , ..., w = w 1 . . . w . Let ϕ c be the ratio between the number of constrains satisfied in the set (14) for a specific solution and its cardinality. Table I shows a comparison between ϕ c and its maximum allowed value for simulations in Fig. 3. The fact that, in many cases, final configurations look like d-configurations might mean a good indication that a global minimum for g is attained in each proposed example, since cost l f o d takes into account the total sum of the potentials for each (i, j) only. On the other hand, it is remarkably worth to recall that this term cannot be generally minimized just by satisfying the maximum number of constraints in (14), as example (2,5) in Fig. 3(a) shows. Indeed, even though a five-point truss is the best shape to achieve, a regular pentagon seems to represent a lower minimum for l f o d , probably because the latter is more compact in terms of inter-agent distances.

C. Complex trajectory tracking
In this subsection, complex desired trajectories are required to be tracked by agents' barycenter while feasible constraints with the same form of (7) are imposed. In 4(a), n = 3 agents follow a curve γ 2 (t) parametrized as γ 2 (t) : where r = 2 m and v = 1 m s −1 . In this framework, robots are also required to gather in a rectangular-triangle formation characterized by desired inter-agent distances equal to (3,4,5) m. Initial conditions for the system are set as follows: Whereas, in Fig. 4(b), n = 4 agents follow a curve γ 3 (t) parametrized as where r = 15 m and v = 2 m s −1 . In this framework, agents are also required to gather in a planar squared formation characterized by desired inter-agent distances equal to 5 m (sides) and 5 √ 2 m (diagonals). Initial conditions for the system are set as follows: In both cases, q p = 10 2 m −2 , q v = 1 m −2 s 2 and M axIter = 50 are assigned, since we have noticed that the value of q p visibly affects the tracking performances 6 , especially in the three-dimensional case.

D. Invariant properties of the solution
One intriguing property appears during this study: initial conditions and desired trajectory could affect the evolution of the system in terms of invariant sets. More precisely, we can state the following Conjecture 4.1: Let S be a subspace of R M such that x B,des (·) ∈ S and x(0) ∈ S. By assuming all the hypotheses stated in Sec. II, the trajectory x(·) computed by Alg. 1 entirely belongs to S. This conjecture is supported by the analysis provided in the forthcoming simulation, where q p = 10 m −2 , q v = 1 m −2 s 2 and M axIter = 50 are assigned and initial conditions for scenarios in Fig. 5(a)

V. CONCLUSIONS AND FUTURE DIRECTIONS
We have formalized a simplified OIFT problem whose minimum-energy solution can be computed by a numerical tool called PRONTO. In order to achieve precise final configurations for the formations we want to drive, a meticulous use of potential functions have been made. In particular, simulations not only show the correctness of this approach but also give a fair description of its robustness when infeasible geometric constraints are imposed to the inter-agent distances. Numerical results also exhibit the versatility of this algorithm, since relatively complex desired trajectories can be tracked by tuning opportunely the weighting parameters contained in a suitable cost function that captures all the requirements.
Possible future directions for this work might be represented by the application of PRONTO to a real OIFT problem with a nonlinear dynamics or the extension to a time-varying formation framework as well as a decentralized paradigm envisaging clusters of mobile robots with relevant communication constraints on the information sharing. Further inspections on the asymptotic stability for a system of agents steered by this optimization approach may likewise depict a fascinating research in this area.