Multicluster Design and Control of Large-Scale Affine Formations

1 Multicluster Design and Control of Lar ge-Scale Af fine F ormations Zhonggang Li, Graduate Student Member , IEEE , Geert Leus, F ellow , IEEE , and Raj Thilak Rajan, Senior Member , IEEE Abstract —Con ventional affine formation control (AFC) em- powers a network of agents with flexible but collective motions - a potential which has not yet been exploited for large-scale swarms. One of the key bottlenecks lies in the design of an interaction graph, characterized by the Laplacian-like str ess matrix. Efficient and scalable design solutions often yield suboptimal solutions on various performance metrics, e.g., conv ergence speed and communication cost, to name a few . The current state-of-the- art algorithms for finding optimal solutions ar e computationally expensive and therefor e not scalable. In this work, we propose a more efficient optimal design for any generic configuration, with the potential to further reduce complexity for a large class of nongeneric r otationally symmetric configurations. Furthermore, we introduce a multicluster control framework that offers an additional scalability improvement, enabling not only collective affine motions as in conventional AFC but also partially inde- pendent motions naturally desired for large-scale swarms. The overall design is compatible with a swarm size of several hundred agents with fast formation conver gence, as compared to up to only a few dozen agents by existing methods. Experimentally , we benchmark the perf ormance of our algorithm compared with several state-of-the-art solutions and demonstrate the capabilities of our pr oposed control strategies. Index T erms —distributed control and optimization, communi- cation networks, networked control systems I . I N T R O D U C T I O N I N recent years, robotic swarms hav e attracted significant attention for their ability to ex ecute complex, large-scale missions [2], [3], [4]. These swarms are deployed across a wide v ariety of domains, from en vironmental exploration [5] and collective payload transport [6], to synchronized satellite constellations [7], some of which rely on tightly coordinated behaviors among the individual robots. For these specific missions, where maintaining a spatial relationship during collectiv e motions is essential, formation control serves as an essential mechanism [8], [9], [10]. A distributed formation control system uses relativ e information such as interagent dis- tances [11], displacements [12], bearings [13], etc., to achiev e and maintain a desired geometric pattern in 2D or 3D space, and possibly engage in continuous formation maneuv ers. Such a system is typically characterized by a framew ork consisting This work was partially presented at the 33rd European Signal Processing Conference (EUSIPCO) in 2025 [1], and is partially funded by the Sensor AI Lab, under the AI Labs program of Delft University of T echnology and by the European Commission Ke y Digital T echnologies Joint Undertaking - Research and Innov ation (HORIZON-KDT -JU-2023-2-RIA), under grant agreement No 101139996, the ShapeFuture project. The authors are with the Signal Processing Systems group, the Faculty of EEMCS, Delft University of T echnology , 2628CD, Delft, The Netherlands (email: { z.li-22, g.j.t.leus, r .t.rajan } @tudelft.nl). of two main components: (a) a graph, where the edges denote the communication links for information exchange, and (b) a configuration that instantiates the graph nodes with spatial coordinates. The design of such a frame work, similar to domains such as sensor networks [14], distrib uted optimization [15], etc., plays a critical role in determining the o verall stability , performance, and communication cost of the system. Recently , affine formation control (AFC) has garnered sig- nificant attention for its maneuverability and coordination flexibility [16], [17], [18]. Unlike traditional rigid forma- tion control, AFC permits continuous geometric transfor- mations—such as scaling, shearing, and rotation—making it adv antageous for adaptive mission planning in dynamic en vironments. Fundamentally , AFC relies on a consensus- based framew ork where control inputs are derived from the relativ e displacements among neighboring agents. While this distributed architecture naturally promotes scalability , existing literature rarely addresses the practical challenges of deploying AFC in large-scale swarms. The core of these challenges is the network topology: instead of a standard graph Laplacian, AFC employs a stress matrix [19], a generalized Laplacian with real-valued edge weights intrinsically tied to graph rigidity . Designing a stress matrix is the foundational step for AFC implementation. While the legitimac y of stabilization is a hard constraint, sev eral additional critical aspects should also be considered for large-scale swarms, such as the sparsity of the underlying graph determining the communication overhead, the con vergence speed of AFC giv en an arbitrary initialization, and the computational complexity for practical feasibility . Designing a stabilizing stress matrix for fast AFC conv er- gence while maintaining a sparse communication topology is challenging, as the communication capacity affects the efficienc y of information exchange and thus influences the con ver gence [20], [21]. Earlier mathematical literature study rigid geometric designs that guarantee the existence of a stress matrix [22], [19]. Subsequent works formulate numerical feasibility problems for computing such matrices [16], [17]. The focus is on finding a feasible stabilizing stress matrix with a prescribed graph topology , rather than optimizing any objecti ve. This prescription can be impractical for large swarms in 2D and 3D. Similar approaches are adopted in [23] and [24], in which a sparsity pattern of the underlying graph is predetermined. Specifically , [23] adopts a sparse nullspace reconstruction approach that sets specific locations in the stress matrix to zero, and calculates the rest to satisfy the constraints. Expanding graphs is also common in constructing rigid frameworks with stress matrices [24], [25], where an 2 initial small base graph is considered, and the other nodes are added sequentially by strategically connecting to a few nodes in the existing graph. Both methods yield substantial sparsification compared with a complete graph, resulting in low communication overhead, and admit analytical solutions that are computationally efficient. Ho wever , they often rely on the genericity assumption of the given node configuration, and the sparsity pattern or the connecting strategy is heuristic. The critical aspect of con vergence speed is not considered either , especially for large networks. A systematic search for an optimal solution with respect to both con vergence and com- munication cost is dev eloped in [26], in which the existence of edges is modeled by binary variables and the problem is formulated as a mixed-inte ger semidefinite program (MISDP). Although this framew ork admits optimized solutions, it raises concerns about the computational complexity and optimality due to the nonconv ex formulation. As such, this motiv ates a stress design method that jointly optimizes the critical objec- tiv es while maintaining manageable computational complexity for large networks. Additionally , the absence of systematic comparisons in the literature calls for a quantitati ve ev aluation of existing methods. In this paper , we present an optimal design of stress matrices through an efficient con vex optimization framework for any giv en configuration. W e focus on network sparsification for efficient communication while guaranteeing AFC con vergence acceleration. Compared with our preliminary results in [1], our additional contributions include the following. • W e reveal the optimal structure of the stress matrix for some special nongeneric geometries, i.e., rotationally symmetric configurations, which further facilitates a re- duction in the computational complexity . • T o ensure the stress design is compatible with lar ger-scale networks, we propose a control framework that partitions the swarm into smaller clusters and aligns their motions through bridging nodes. • W e present insights on the stability and provide condi- tions when the clusters exhibit collectiv e behaviors, as well as when inter-cluster flexibility is introduced. • W e perform a comprehensi ve performance benchmark of our proposed stress design solution against sev eral state- of-the-art solutions across generic and nongeneric test cases ranging from small-scale ( ≤ 10 nodes) to large- scale ( ≥ 100 nodes) setups in 2D and 3D. The rest of the paper is organized as follows. In Section II, we introduce the preliminaries of graph rigidity and stress- based AFC. In Section III, we propose a con ve x optimization framew ork for the stress design and discuss the choice of hyperparameters to balance the objecti ves. This is considered a generic solution compatible with any given configuration. For symmetric configurations, we discuss their optimal structure under our proposed formulation, and propose the complexity- reduced unique stress identifier (USI) method in Section IV. In Section V, we introduce the multicluster control framew ork based on a classical control law , discuss its stability , and provide conditions under which the clusters exhibit fully or partially collecti ve beha viors. Finally , we compare our pro- posed solutions with the state-of-the-art through a comprehen- siv e benchmark and demonstrate the capability of multicluster control in Section VI. Concluding remarks and future work are briefly outlined in Section VII. Notations. V ectors and matrices are represented by lower- case and uppercase boldface letters, respectiv ely , such as a and A . Sets and graphs are represented using calligraphic letters, e.g., A . V ectors of length N of all ones and zeros are denoted by 1 N and 0 N , respecti vely , and their matrix versions are similarly 1 M × N and 0 M × N . An identity matrix of size N is denoted by I N . The Kronecker product is ⊗ and the vectorization operator is vec( · ) . The diag ( · ) operator creates a diagonal matrix from a vector . Moreov er , tr( · ) denotes the trace operator and λ k ( A ) denotes the k -th smallest eigen v alue of a symmetric matrix A . W e use [ A ] ij and [ a ] i to denote the elements of matrices and vectors, respectively . Units are enclosed by square brackets, e.g., seconds [s]. I I . F U N DA M E N TA L S A. Graphs and Rigidity Theory Consider N mobile agents in D -dimensional Euclidean space where N ≥ D + 1 . They are interconnected with a communication network, modeled by an undirected graph G = ( V , E ) , where the vertices V = { 1 , ..., N } denote the agents, and the edges E ⊆ V × V denote the communication links. W e use N = |V | and M = |E | as shorthand for the number of ver- tices and edges, respectiv ely . The set of neighbors of a node i is defined as N i = { j ∈ V : ( i, j ) ∈ E } . Each node i ∈ V is as- signed a position p i ∈ R D , and the configuration of the nodes in V is P ( V ) = [ p 1 , . . . , p N ] ∈ R D × N . W e also define an augmented configuration ¯ P ( V ) =  P ⊤ , 1 N  ⊤ ∈ R ( D +1) × N . They will be shorthanded to P and ¯ P respecti vely , unless intended for particular sets of vertices. A g eneric configuration [27] has algebraically independent node coordinates, i.e., there are no geometric constraints among nodes. A configuration of random node coordinates has probability one of being generic. T ypical nongeneric configurations include grids, a regular polygon, etc. A frame work F = ( G , P ) pairs the graph G with its configuration P . T wo frameworks ( G , P ) and ( G , P ′ ) are equivalent if ∥ p i − p j ∥ 2 =   p ′ i − p ′ j   2 for all ( i, j ) ∈ E , and congruent if this distance equality holds for all node pairs ( i, j ) ∈ V × V . A frame work F is globally rigid (often termed rigid ) if ev ery equiv alent framework is also congruent. Furthermore, F is universally rigid if it remains globally rigid in any higher-dimensional embedding R D ′ where D ′ ≥ D . Fig. 1 illustrates these rigidity concepts. Univ ersal rigidity can be algebraically represented by an equilibrium stress ω ij ∈ R for ev ery edge ( i, j ) ∈ E , which together satisfy X ( i,j ) ∈E ω ij ( p i − p j ) = 0 D . (1) A more compact form of (1) is Ω P ⊤ = 0 N × D , where Ω ∈ R N × N is the str ess matrix defined as [ Ω ] ij =      0 , if ( i, j ) / ∈ E − ω ij , if i  = j, ( i, j ) ∈ E P k ∈N i ω ik if i = j . (2) 3 (a) flexible (b) globally rigid (c) universally rigid Fig. 1. Examples of frameworks in R 2 with increasing rigidity [16] Alternativ ely , the equilibrium stress and the stress matrix are linked in matrix-v ector form through the graph incidence matrix B ∈ R N × M Ω = B diag( ω ) B ⊤ , (3) where ω ∈ R M is a vector containing all the equilibrium stresses. Note that Ω reduces to a standard graph Laplacian if diag( ω ) = I M , i.e., equal weights for the edges, and it also has 1 N in the nullspace like the Laplacian. The following theorem establishes the important properties of Ω related to univ ersal rigidity . Theorem 1. ( Universally Rigid F rameworks and Str ess Matri- ces ) Giv en a frame work F = ( G , P ) with P being a generic configuration, F is univ ersally rigid if and only if there e xists a stabilizable stress matrix Ω , i.e., a positiv e semidefinite (PSD) stress matrix with rank N − D − 1 . Proof . See [28], [17], [26]. ■ Note that a generic configuration is typically assumed by the literature for mathematical guarantees, b ut nongeneric geome- tries are also useful for formation control applications. It is worth mentioning that our proposed stress design approaches work for both generic and nongeneric configurations. B. Single-Cluster Affine F ormation Contr ol Giv en a stress matrix Ω under a universally rigid frame work, affine formation control (AFC) tries to steer the real-time position z i ∈ R D into a target position p i for agent i ∈ V . T o achie ve this, we assume agents adopt single-inte grator dynamics ˙ z i = u i where u i is a velocity control input. The con ventional consensus-based control law is designed as [16] u i = − X j ∈N i [ Ω ] ij ( z j − z i ) , (4) which has been extended and enriched for se veral scenarios [17], [29]. The collecti ve dynamics of (4) can be described by the linear time-in variant (L TI) system ˙ z = − ( Ω ⊗ I D ) z , (5) where z =  z ⊤ 1 , ..., z ⊤ N  ⊤ ∈ R DN . Recall that the stress matrix from (1) satisfies Ω P ⊤ = 0 , which can be re written as ( Ω ⊗ I D ) p = 0 , where p = v ec( P ) . It is then straightforward to see that p is an equilibrium of system (5). This is considered an extension of the av erage consensus system using the graph Laplacian, which has only a 1 -dimensional nullspace corresponding to the 1 N vector , i.e., only configurations with all nodes in the same location can be an equilibrium. Moreov er , since Ω is positive semidefinite, system (5) is globally and exponentially stable [17], meaning that the formation asymptotically conv erges to the equilibrium space containing p given any initialization. Note that the targeted equilibrium p can be reached up to an affine transformation (hence its usefulness for AFC), and can be uniquely deter- mined giv en a minimum of D + 1 anchor nodes, which are commonly named leaders. Conceptually , they tie up the loose degrees of freedom from the flexible affine transformations and thus can guide continuous formation maneuvers. This is fa vored since only a small set of leaders needs to be manipu- lated, while the majority remains fully autonomous. Ho wever , the collectiv e affine flexibility may still be considered limited in some cases for large-scale formation as a whole. C. Pr oblem F ormulation Giv en a potentially large set of agents V ( N ≥ 100 ) with any configuration P , we primarily seek a stabilizable stress matrix that yields fast con ver gence for system (5) and a sparse structure for lower communication load. In consensus theory , the second smallest eigen value of the Laplacian (Fiedler value) governs con vergence speed [21], as it determines the slowest decaying component of the error . Prior works hav e aimed to maximize this eigen value for faster conv ergence [20]. Similarly , we aim to maximize the smallest nonzero eigen value, i.e., the ( D + 2 )-th smallest eigenv alue of Ω or λ D +2 ( Ω ) for stabilizable stress matrices. Additionally , we seek a control framew ork that preserves the affine flexibility but could introduce additional degrees of freedom beyond collectiv e affine transformations. I I I . O P T I M A L T O P O L O G Y D E S I G N In this section, we focus on optimal stress matrix design using a con vex formulation, considering all the nodes as one group. W e first present a preliminary formulation that directly reflects the objecti ves and constraints b ut is intractable. W e then reformulate it into a tractable con vex optimization prob- lem and provide insights into the introduced hyperparameters. Recall from (3) that the stress matrix can be constructed using the graph incidence matrix giv en a stress vector . W e initialize the graph as a complete graph with incidence matrix ¯ B ∈ R N × ¯ M having ¯ M = N ( N − 1) 2 edges. W e then aim to find a sparse vector ¯ ω ∈ R ¯ M such that Ω is sparse by (3). The effecti ve number of edges will be M = ∥ ¯ ω ∥ 0 . Hence we pose the following problem P 0 : P 0 : minimize ¯ ω , Ω ∥ ¯ ω ∥ 0 − α λ D +2 ( Ω ) (6a) subject to Ω = ¯ B diag ( ¯ ω ) ¯ B ⊤ (6b) rank( Ω ) = N − D − 1 (6c) Ω ⪰ 0 (6d) Ω ¯ P ⊤ = 0 (6e) where α is a weighting parameter between the network sparsity by ∥ ¯ ω ∥ 0 and the conv ergence speed by λ D +2 ( Ω ) . As can be observed, P 0 is difficult to solve due to the L0-norm and the 4 rank constraint. In the next section, we con vexify this problem to a tractable semidefinite program (SDP). A. Con vexification and Eigen value Maximization A common con vex relaxation for promoting sparsity re- places the L0-norm with an L1-norm, i.e., we relax ∥ ¯ ω ∥ 0 to ∥ ¯ ω ∥ 1 in the objective. Further, assuming a full-rank ¯ P , we let Q ∈ R N × ( N − D − 1) with orthonormal columns span the kernel of ¯ P . Then, (6c), (6d), and (6e) imply that Q ⊤ Ω Q ∈ R ( N − D − 1) × ( N − D − 1) is a positiv e definite (PD) matrix, since the null-subspace of Ω is projected out by Q , i.e., Q ⊤ Ω Q ≻ 0 . As such, this PD constraint is equiv alent to (6c) under (6d) and (6e). T o maximize λ D +2 ( Ω ) , we explicitly use the trace of the stress matrix tr( Ω ) = tr  Q ⊤ Ω Q  , a con vex proxy which effecti vely sums all the (nonzero) eigen values of Ω , together with a spectral norm constraint ∥ Ω ∥ 2 ≤ β , where we limit the largest eigen value of Ω to β > 0 . W ithout additional objectiv es and constraints, λ D +2 ( Ω ) is ef fectiv ely maximized to λ D +2 ( Ω ) = ... = λ max ( Ω ) = β with a minimum condition number κ ( Ω ) = 1 . In our formulation, this boosting effect is further shaped by the sparsity term and other constraints. It is worth mentioning that trace regularization is also commonly seen in network design literature, such as [20]. Additionally , a low condition number has been shown to improve the robustness of AFC against time delays [26] and hence is aligned with our objectiv e for conv ergence acceleration. Having the above remodeling, a working con vex formula- tion is shown in P 1 P 1 : minimize ¯ ω , Ω f ( ¯ ω ) = ∥ ¯ ω ∥ 1 − α tr( Ω ) (7a) subject to Ω = ¯ B diag ( ¯ ω ) ¯ B ⊤ (7b) ∥ Ω ∥ 2 ≤ β (7c) Q ⊤ Ω Q − γ I N − D − 1 ⪰ 0 (7d) Ω ¯ P ⊤ = 0 (7e) where α > 0 is the weighting parameter of the objectives and β > 0 upper bounds the spectral norm of Ω . W e also set a lower bound γ > 0 for (7d) for two explicit reasons, (a) positiv e-definiteness Q ⊤ Ω Q ≻ 0 cannot be strictly enforced in numerical solvers; and (b) when the L1 term is dominating, we want to av oid the trivial solution ¯ ω = 0 . Note that P 1 can be solved by off-the-shelf solvers, but could be cumbersome because the problem scales badly , as both ¯ ω and Ω are optimization variables while linked through an equality constraint. Next, we further simplify P 1 by using only ¯ ω as the optimization variable, yielding an equiv alent but efficient final formulation. B. Pr oposed F ramework Substituting (7b) into (7e), we hav e ¯ P ¯ B diag ( ¯ ω ) ¯ B ⊤ = ¯ P ¯ B diag ( ¯ ω )  ¯ b 1 , ..., ¯ b N  = 0 where ¯ b i ∈ R ¯ M , ∀ i ∈ V is the i -th column of ¯ B ⊤ . Observe that ¯ P ¯ B diag ( ¯ ω ) ¯ b i = ¯ P ¯ B diag  ¯ b i  ¯ ω = 0 . As such, we can construct a matrix E ∈ R N ( D +1) × ¯ M with structure E =    ¯ P ¯ B diag  ¯ b 1  . . . ¯ P ¯ B diag  ¯ b N     , (8) such that E ¯ ω = 0 , which can replace constraint (7e). Ne xt, we combine (7b) and (7d) such that Q ⊤ Ω Q = Ψ diag ( ¯ ω ) Ψ ⊤ ≻ 0 , where Ψ = Q ⊤ ¯ B . W e can then also write tr  Q ⊤ Ω Q  = tr  Ψ diag ( ¯ ω ) Ψ ⊤  = tr  Ψ ⊤ Ψ diag ( ¯ ω )  (9) = tr  diag  Ψ ⊤ Ψ  ¯ ω  = ψ ⊤ ¯ ω , where ψ = diag  Ψ ⊤ Ψ  . As such, we could re write P 1 into an equiv alent form P 2 P 2 : minimize ¯ ω ∥ ¯ ω ∥ 1 − α ψ ⊤ ¯ ω (10a) subject to Ψ diag ( ¯ ω ) Ψ ⊤ − γ I N − D − 1 ⪰ 0 (10b)   ¯ B diag ( ¯ ω ) ¯ B ⊤   2 ≤ β (10c) E ¯ ω = 0 (10d) where there is only ¯ ω as optimization variable. Next, we discuss the choice of hyperparameters, namely α , β , and γ in the formulation. C. Choice of Hyperparameters The need for β > γ is straightforward to ensure feasibility since the y represent the largest and smallest eigen values of Ω , respectively . They should also be sufficiently apart for a relativ ely large feasible re gion. The objectiv e function contains a sparsity term ∥ ¯ ω ∥ 1 (geometrically a cone) and a con ver gence term − α ψ ⊤ ¯ ω (a hyperplane). If the con ver gence term is dominant, i.e., α is large, the hyperplane will unfold the cone such that its global minimum vanishes, in which case the solution is minimized at the boundary of the feasible region, which is not necessarily a sparse solution. Mathematically , to ensure the existence of a global minimum of the objective function, 0 ¯ M should be in the subdif ferential of the objective function, i.e., 0 ¯ M ∈ ∂  ∥ ¯ ω ∥ 1 − α ψ ⊤ ¯ ω  , (11) which directly yields that all elements of α ψ are in the range ( − 1 , 1) . Recalling that all elements of ψ are nonnegati ve from the definition ψ = diag  Ψ ⊤ Ψ  and that α > 0 , we could thus define a critical v alue α ∗ = 1 / ∥ ψ ∥ ∞ , where the infinity norm is the largest value of ψ . The solution will be sparse with maximally accelerated con vergence at the critical value α ∗ , and the sparsity will trade of f with faster con ver gence for larger α . Note that values below α ∗ do not necessarily yield a sparser solution. Rather , the dominant L1-norm forces all edge weights closer to zero, including the activ e stress values, which mak es it dif ficult to e xtract the graph topology by setting a threshold. 5 (a) An octagon in 2D (b) An octahedron in 3D. (c) A cuboctahedron in 3D. Fig. 2. A few examples of rotationally symmetric configurations. Note that the optimal stress for these configurations also applies to their affine transformations. I V . T H E U N I Q U E S T R E S S I D E N T I FI E R ( U S I ) M E T H O D The previous section proposes a con ve x formulation to design the stress matrix for any given configuration. In this section, we narrow our scope down to a class of highly symmetric configurations in 2D and 3D, namely , rotationally symmetric configurations, with se veral e xamples sho wn in Fig. 2. W e rev eal that the optimal stress matrices for these configurations are highly structured, and the stress v alues are also highly repetiti ve. Therefore, the number of distinct values in the stress for optimization significantly decreases, which motiv ates the design of a shortcut formulation that optimizes only the unique values of the stress with lower computational complexity . W e name this approach the unique stress identifier (USI) method. A. Optimal Str ess for Symmetric Configurations Rotationally symmetric configurations are tightly coupled with node permutations, which is a crucial property we lev er- age to re veal the optimal stress structure. These configurations are formalized in the following definition. Definition 1 ( Rotationally symmetric configurations ) . A rota- tionally symmetric configuration P ∈ R D × N satisfies R k P = P Π k , (12) where R k ∈ SO ( D ) and Π k ∈ { 0 , 1 } N × N , for k = 1 , ..., K , is the k -th rotation-permutation pair that preserves the config- uration. This means that rotationally symmetric configurations are in variant to rotations up to a reordering of the vertices. Since the stress matrix Ω coexists with the stress vector ω through the graph incidence matrix (3), there exists a corresponding edge permutation ¯ Π E ,k ∈ { 0 , 1 } ¯ M × ¯ M for the complete graph characterized by the incidence matrix ¯ B to a giv en node permutation Π k , such that Π k ΩΠ k ⊤ = ¯ B diag  ¯ Π E ,k ¯ ω  ¯ B ⊤ . (13) This is considered the equiv alence of permutations in the node and edge domain. Next, we build up towards the structure of optimal stress matrices under symmetric configurations. For mathematical simplicity , we focus on P 1 (7) instead of the equiv alent P 2 (10). Definition 2 ( P ermutation in variance of feasible sets of P 1 ) . Denote the feasible set of P 1 as C . Giv en a rotationally symmetric P as in Definition 1, then the feasible set C is permutation in variant if, for any feasible ¯ ω ∈ C , all the K permutations are also feasible, i.e., ¯ Π E ,k ¯ w ∈ C , k = 1 , ..., K , where ¯ Π E ,k ∈ { 0 , 1 } ¯ M × ¯ M is an edge permutation related to the k -th symmetry rotation. Definition 3 ( P ermutation in variance of objectives of P 1 ) . De- fine the objective function of P 1 as f ( ¯ ω ) , where f : R ¯ M 7→ R . Giv en a rotationally P in Definition 1, f ( ¯ ω ) is permutation in variant if the permutation of the argument ¯ ω results in the same objectiv e value, i.e., f ( ¯ Π E ,k ¯ ω ) = f ( ¯ ω ) , k = 1 , ..., K . Lemma 1 ( P ermutation in variance of P 1 ) . The objective function (7a) and the feasible set C defined by constraints (7b) to (7e) of P 1 are both permutation in variant. Proof . See Appendix A. ■ Lemma 1 establishes that if a stress matrix for rotationally symmetric configurations is feasible under P 1 , then applying permutations to the stress will preserve the feasibility and the cost value. W e then define the permutation-in variant stress matrices, which we show are optimal under P 1 . Definition 4 ( P ermutation-invariant str ess matrices Ω pi ) . For a given rotationally symmetric configuration, a permutation- in variant Ω pi ∈ R N × N with its vector counterpart ¯ ω pi ∈ R ¯ M satisfying Ω pi = ¯ B diag ( ¯ ω pi ) ¯ B ⊤ , satisfies Ω pi = Π k Ω pi Π k ⊤ , (14) where Π k for k = 1 , ..., K denotes the permutation from (12) in Definition 1. Theorem 2 ( Optimality of permutation-in variant stress ) . Giv en a rotationally symmetric configuration, there exists a permutation-in variant Ω pi that is feasible and optimal for P 1 . Proof . Assume ¯ w ∗ ∈ R ¯ M is a minimizer of P 1 with its matrix counterpart Ω ∗ = ¯ B diag ( ¯ ω ∗ ) ¯ B ⊤ , which is not necessarily permutation-in variant. Using the in variance of feasible sets and objective functions in Lemma 1, we can find K feasible stress by permuting the configuration with the same objective value, i.e., ¯ Π E ,k ¯ ω ∗ or Π k Ω ∗ Π ⊤ k for k = 1 , ..., K . Then we construct a permutation-inv ariant stress by av eraging across all the permutations, i.e., Ω pi = 1 K P K k =1 Π k Ω ∗ Π k ⊤ and ¯ ω pi = 1 K K X k =1 ¯ Π E ,k ¯ ω ∗ , (15) which is also feasible since it is a con vex combination of feasible points. Using Jensen’ s inequality for con vex functions, observe that f ( ¯ ω pi ) = f 1 K K X k =1 ¯ Π E ,k ¯ ω ∗ ! ≤ 1 K K X k =1 f  ¯ Π E ,k ¯ ω ∗  . (16) Since the permutations have the same objectiv e from the in variance of P 1 , we know f  ¯ Π E ,k ¯ ω ∗  = f ( ¯ ω ∗ ) for all k , which can be substituted into the inequality and yield f ( ¯ ω pi ) ≤ 1 K K X k =1 f ( ¯ ω ∗ ) = 1 K K f ( ¯ ω ∗ ) = f ( ¯ ω ∗ ) . (17) This means a permutation-in variant stress will achie ve equal or lower cost than the assumed minimizer with arbitrary structure. 6 In other words, the permutation-in variant stress is the optimal solution to P 1 . ■ From Theorem 2 it is guaranteed that problem P 1 (or equiv- alently P 2 ) under a rotationally symmetric configuration P following (12), has an optimal permutation-in variant solution Ω pi . Next, we show that Ω pi has identical stress values if their corresponding edges in the complete graph can be mapped from one to the other using a symmetric rotation-permutation pair , just like the in variance (14) suggests. Corollary 1 ( Edge equivalence in the stress matrix ) . Let P be a rotationally symmetric configuration and Ω pi satisfying (14) be the optimal solution. T wo verte x pairs, ( i, j ) and ( g , l ) , are defined as symmetrically equivalent if there exists a symmetry rotation-permutation pair ( R k , Π k ) that maps vertex i to g and verte x j to l . If the vertex pair ( i, j ) is symmetrically equiv alent to the verte x pair ( g , l ) , their corresponding stresses are equal, i.e., [ Ω pi ] ij = [ Ω pi ] g l . (18) Proof . If permutation Π k maps index i to g and j to l , then [ Π k ] g i = 1 and [ Π k ] lj = 1 . Since Π k is a symmetric permutation, it holds that Ω pi = Π k Ω pi Π k ⊤ , which can be written in its matrix element form as  Π k Ω pi Π ⊤ k  g l = N X p =1 N X q =1 [ Π k ] g p [ Ω pi ] pq  Π ⊤ k  q l . (19) Since [ Π k ] g p = 1 only when p = i and otherwise is zero, and  Π ⊤ k  q l = [ Π k ] lq = 1 only when q = j and otherwise is zero, (19) simplifies to  Π k Ω pi Π ⊤ k  g l = [ Ω pi ] ij = [ Ω pi ] g l . Hence (18) is proved. ■ Corollary 1shows that for symmetric configurations, edges that map from one to the other have the same value in the optimal solution, which motiv ates us to use the distinct edge (stress) classes instead of for the complete graph for the optimization. Before we implement this observ ation into an efficient algorithm design, we discuss a special case, where the nodes lie on a circle in 2D, forming regular polygons. Remark 1 ( Stress structure for cir cular configurations ) . If the nodes lie on a circle and their configuration P is ordered sequentially , the permutation matrices Π k for k = 1 , ..., K for (12) are limited to cyclic permutations. Then the Ω pi is a circulant (also T oeplitz) matrix, as shown in the example in Fig. 3. B. Efficient Algorithm Design Motiv ated by Corollary 1, we would first identify S unique stress classes from the symmetric configuration. Howe ver , trav ersing the matching edge pairs is nontrivial for complicated configurations. It is straightforward to see that having identical edge lengths (Euclidean distances between vertices) is a nec- essary condition for edge symmetric equiv alence. As such, we propose an approximation using the Euclidean distance matrix (EDM) defined by [edm( P )] ij = ∥ p i − p j ∥ 2 2 i ∈ V , j ∈ V , (20) as an initial attempt to identify the S unique stress classes. Edges with the same length shall be considered the same class. After identifying the unique edge classes, we could create a mapping ¯ ω = S ω where ω ∈ R S is the reduced stress vector , and S ∈ { 0 , 1 } ¯ M × S is a selection matrix with s s ∈ { 0 , 1 } ¯ M denoting the s -th column. Then the mapping could be expanded to ¯ ω = P S s =1 ω s s s where ω s is the s-th element of ω . The objective of P 2 (10a) could be substituted into ∥ ¯ ω ∥ 1 − α ψ ⊤ ¯ ω = ∥ S ω ∥ 1 − α ψ ⊤ ( S ω ) (21a) = S X s =1 ( 1 ⊤ ¯ M s s ) ω s − α ( S ⊤ ψ ) ⊤ ω (21b) = c ⊤ | ω | − α ψ ⊤ ω , (21c) where c = S ⊤ 1 ¯ M ∈ Z S + and ψ = S ⊤ ψ ∈ R S can be precomputed. The PSD constraint (10b) can be replaced by Ψ diag ( ¯ ω ) Ψ ⊤ = Ψ diag S X s =1 ω s s s ! Ψ ⊤ (22a) = S X s =1 ω s  Ψ diag ( s s ) Ψ ⊤  (22b) = S X s =1 ω s Ψ s , (22c) where Ψ s = Ψ diag ( s s ) Ψ ⊤ can also be precomputed. Simi- larly , the spectral norm constraint (10c) is   ¯ B diag ( ¯ ω ) ¯ B ⊤   2 = S X s =1 ω s  ¯ B diag ( s s ) ¯ B ⊤  (23a) = S X s =1 ω s M s , (23b) where the basis matrix M s = ¯ B diag ( s s ) ¯ B ⊤ is kno wn. Lastly , the nullspace constraint (10d) can be expressed as E ¯ ω = ( E S ) ω = E ω = 0 , (24) where E = E S . Having reduced all the expressions in P 2 , we introduce the reduced formulation P 3 P 3 : minimize ω ∈ R S c ⊤ | ω | − α ψ ⊤ ω (25a) subject to S X s =1 ω s Ψ s − γ I N − D − 1 ⪰ 0 (25b)      S X s =1 ω s M s      2 ≤ β (25c) E ω = 0 (25d) A clarifying diagram is provided in (26), summarizing all the proposed formulations for the stress design. P 0 relax = = = ⇒ P 1 equiv alent ⇐ = = = = = ⇒ P 2 symmetric ⇐ = = = = = = = ⇒ configurations P 3 (26) Fig. 3 visualizes a few e xamples of the stress matrix optimized through P 3 and P 2 in comparison with their corresponding EDMs. As can be seen, the EDM correctly identifies the stress 7 Generic USI EDM Fig. 3. Stress matrices compared with EDMs using an octagon in Fig. 2(a) (upper) and a cuboctahedron in Fig. 2(c) (lo wer), which illustrate the equiv alence of the proposed optimization (Generic P 2 ) and the unique stress identifier (USI) method P 3 on rotationally symmetric configurations. The stress matrices for the octagon are symmetric, T oeplitz, and circulant. The diagonals of matrices are made empty as they do not represent distinct edge classes but are determined by the others. Matching colors represent the same stress values for the first two columns, and the same stress classes to compare with the EDM. classes in comparison with the generic formulation P 2 , and the reduced-size formulation P 3 can find identical solutions. In practice, the EDM approximation often works robustly for symmetric configurations, although it identifies equal or fewer stress classes than the rigorous equiv alence requirement. C. Complexity Analysis In general, semidefinite programming can be solv ed by interior point methods (IPMs) in roughly cubic complexity w .r .t. the number of decision variables [30], i.e., the size of the stress vector ¯ ω for P 2 , or ω for P 3 . This is already advantageous to the NP-hard MISDP formulation [26] that scales exponentially . W e give a complexity analysis of the proposed P 2 and P 3 w .r .t. the number of decision variables and the prospects of further reduction using the multicluster framew ork in the following sections. Recollect that the generic formulation P 2 has ¯ M = N ( N − 1) 2 variables for N nodes, which scales quadratically O ( N 2 ) in terms of the number of variables. For symmetric configura- tions, this reduces to O ( S ) using P 3 . For regular 2D polygons, S scales linearly with the number of nodes, hence admitting O ( N ) v ariables. For 3D rotationally symmetric structures, S should be discussed on a case-by-case basis. Howe ver , for highly symmetric configurations, such as the Archimedean solids, it is typical that S ≪ N 2 , which permits a massiv e reduction in the complexity . An e xample is the truncated icosahedron with 60 nodes (used in our benchmark as sho wn in Fig. 6 (d)), which has ¯ M = 1770 v ariables for P 2 but has only S = 21 variables using our USI for P 3 . Although the rotationally symmetric configurations and their af fine transformations cover a large class of nongeneric geometries, P 2 for generic configurations is still inevitable in man y cases. The quadratic scaling in SDP variables can be concerning for a large number of nodes. The next section introduces the multicluster control framew ork, where the large number of nodes is partitioned into smaller clusters, whose stress matrices can be designed independently and more ef fi- ciently , which further reduces the overall design complexity . V . M U LT I C L U S T E R C O N T RO L O F L A R G E N E T W O R K S In previous sections, we established a con vex optimization framew ork for designing sparse networks with high con ver- gence speed. The generic formulation P 2 scales quadratically with the nodes for an SDP formulation, which becomes de- manding for large-scale networks ( N > 100 ). In this section, we introduce a di vide-and-conquer strategy that partitions a large network into sev eral smaller subnetworks with ov erlap- ping nodes, thereby significantly reducing the ov erall design complexity . Additionally , long-distance communications can also be avoided with the overlapping nodes acting as relays. W e propose control strate gies with proofs of stability and analyze the associated conv ergence properties. Furthermore, we deri ve conditions on the o verlapping nodes under which the entire ensemble executes a unified collective af fine motion like the single cluster case, as well as conditions under which dif- ferent clusters exhibit coordinated collectiv e behaviors while preserving individual flexibility . W e emphasize that our pro- posed framework is instantiated using the classic control law (4), but it also applies to its generalizations. Formally , given a nominal configuration P ( V ) ∈ R D × N , we divide it into C subconfigurations P c ( V c ) ∈ R D × N c for c = 1 , ..., C , where V c is the subset of nodes for subcon- figuration c , and N c ≥ D + 1 is the number of nodes in subconfiguration c . W e require overlapping nodes among the subconfigurations, formally named bridging nodes , such that P C c =1 N c > N . For each subconfiguration P c , a stress matrix Ω c ∈ R N c × N c can be designed using the generic formulation P 2 . Then the total number of variables is reduced from O ( N 2 ) to O ( P C c =1 N 2 c ) = O ( C N 2 c ) if all clusters are assigned an equal number of nodes N c for a sequential design procedure. If stresses for clusters are designed in parallel, this further reduces to O (max c N 2 c ) . If P c is symmetric, we could use the accelerated P 3 . Let C i denote the number of clusters node i belongs to, then C i = 1 if node i is a nonbridging node, whereas 1 < C i ≤ C if it is a bridging node. Based on these clusters, we can now propose a distributed control strategy as detailed in Algorithm 1. For agent i ∈ V at time t , randomly select a cluster c among all possible C i clusters. Compute the velocity control input for node i related to the selected cluster as u i = − 1 π c i X j ∈N c i [ Ω c ] ij ( z j − z i ) , (27) where N c i is the neighbor set for agent i in cluster c , and the gain π c i is the probability of agent i choosing cluster c . W e opt for a discrete uniform distrib ution, i.e., π c i = 1 C i . For the non- bridging nodes, the selection is deterministic since C i = 1 . A. Stability Analysis Controller (27) represents a linear time-varying (L TV) sys- tem, for which the instantaneous system matrix is time-v arying and asymmetric. W e consider the mean dynamics E [ ˙ z i ] = E [ u i ] . (28) 8 Algorithm 1 Multicluster contr ol 1: Input : Configuration P ∈ R D × N 2: Initialize : z i for i ∈ V 3: Divide P into P c ∈ R D × N c for c = 1 , ..., C 4: Calculate Ω c ∈ R N c × N c for c = 1 , ..., C 5: for i ∈ V do 6: Randomly choose cluster c 7: Exchange for z j ∈ N c j 8: Calculate u i = − 1 π c i P j ∈N c i [ Ω c ] ij ( z j − z i ) 9: Execute ˙ z i = u i 10: end for T aking the expectation for the local dynamics (27), we obtain E [ u i ] = C X c =1 π c i   − 1 π c i X j ∈N c i [ Ω c ] ij ( z j − z i )   (29a) = − C X c =1 X j ∈N c i [ Ω c ] ij ( z j − z i ) . (29b) Notice that the probability term π c i cancels with the gain compensation factor . Re writing this in global form, the mean dynamics become: E [ ˙ z ] = − ( C X c =1 ˜ Ω c ) ⊗ I D ! E [ z ] = − ( Ω ⊗ I D ) E [ z ] (30) where ˜ Ω c ∈ R N × N is the expansion of the stress matrix Ω c into the global dimension with zero padding, and the equiv alent global ensemble Ω = P C c =1 ˜ Ω c . A visualization of (30) is shown in Fig. 4, where the 2 clusters are color coded, and the layers distinguish the stress matrices to be added for the ensemble Ω . Note that the mean dynamics yield a linear time-in variant (L TI) system, which provides insight for the con ver gence property . Theorem 3 ( Con verg ence of mean dynamics ) . Gi ven p = v ec( P ) as a target configuration, the mean dynamics (30) exponentially conv erge to a space containing p given arbitrary initialization. Proof . W e first sho w that the nominal configuration p is still an equilibrium for the equiv alent ensemble stress matrix Ω . The stress matrices Ω c for cluster c = 1 , ..., C satisfy ( Ω c ⊗ I D ) vec( P c ) = 0 . It also satisfies that ( ˜ Ω c ⊗ I D ) p = 0 due to the zero-padding. Then the ensemble Ω , the sum of ˜ Ω c ov er c = 1 , ..., C satisfies ( Ω ⊗ I D ) p = 0 . Now ˜ Ω is also a positiv e-semidefinite (PSD) matrix since it is the sum of C PSD matrices ˜ Ω c . As a result, the dynamics (30) exhibit exponential con vergence. ■ Theorem 3 establishes the stability of the mean dynamics, meaning a non-increasing mean error . W e further giv e insights into the instantaneous dynamics. Remark 2 ( Con ver gence of instantaneous dynamics ) . Let σ ( t ) denote an arbitrary switching pattern representing the agents’ random cluster selections at time t , and let H σ ( t ) ∈ R N × N be the corresponding instantaneous global system matrix, which is a row-wise assembly of ˜ Ω c . The continuous-time system dynamics are ˙ z ( t ) = −  H σ ( t ) ⊗ I D  z ( t ) . (31) Note that H σ ( t ) is not necessarily a PSD matrix to guarantee a nonincreasing error energy using a common L yapunov function V ( z ) = 1 2 ∥ z ( t ) − z ∗ ∥ 2 2 where z ∗ ∈ A ( P ) is a target configuration in the af fine image A ( P ) . This means that at certain time instants, the agents may temporarily pull each other in conflicting directions, thereby increasing the overall error energy , although this has not been experimentally sho wn to be an issue. Note that Theorem 3 only establishes the con ver gence into the space containing the giv en configuration p . Ho wever , whether this space is only spanned by the af fine transformation of p needs further in vestigation, as it determines whether the clusters behav e as one affine maneuvering unity . B. Conditions for Collective Motions In this section, we giv e some insights into the conditions on the bridging nodes for different clusters to adopt collecti ve affine motions, i.e., local affine transformations of each cluster are identical and are thus global. A direct criterion is to ev aluate if the rank of the ensemble stress matrix Ω in (30) is equal to N − D − 1 . A lower rank implies the introduction of additional flex es beyond af fine motions. Howe ver , sho wing the rank of the Ω illustrated in Fig. 4 is not straightforward. Therefore, we show an alternative view of the problem. Theorem 4 ( Conditions on the bridging nodes for collective motions ) . Suppose node i ∈ V c admits z i = Θ c p i + t c =  Θ c t c   p i 1  , (32) where Θ c ∈ R N c × N c and t c ∈ R N c describe the local affine transformations of cluster c . Define the set of bridging nodes B cc ′ = V c ∩ V c ′ for any two connected clusters c, c ′ = 1 , ..., C . Then local affine transformations are identical, i.e., Θ c = Θ c ′ and t c = t c ′ , if the augmented configuration of the bridging nodes ¯ P ( B cc ′ ) ∈ R ( D +1) ×|B cc ′ | is full row-rank. Proof . T o hav e identical local transformations requires z i =  Θ c t c   p i 1  =  Θ c ′ t c ′   p i 1  , ∀ i ∈ B cc ′ , (33) which yields  Θ c t c  −  Θ c ′ t c ′   p i 1  = 0 D , ∀ i ∈ B cc ′ . (34) Stacking for all i ∈ B cc ′ giv es ∆ cc ′ ¯ P ( B cc ′ ) = 0 ( D +1) ×|B cc ′ | , where ∆ cc ′ =  Θ c t c  −  Θ c ′ t c ′  . Then if ¯ P ( B cc ′ ) is full row-rank, we get the trivial solution ∆ cc ′ = 0 . ■ Geometrically , affine motions for two clusters c and c ′ are always aligned if there are at least D + 1 bridging nodes that fully span D -dimensional space for each cluster . If leaders are implemented to guide affine maneuvers, they can be added as required by [17] with no additional restrictions on the distribution to the clusters. 9 cluster 1 cluster 2 bridging nodes E [ ˙ z ] = ! C c =1 ˜ Ω c E [ z ] Ω 1 Ω 2 zero padding Fig. 4. An illustration example of the mean dynamics model (30) with 6 nodes ov er 2 clusters. The dimension is set to D = 1 for simplicity . C. Con ver gence Speed of the Ensemble Graph Since the ensemble graph works as a single cluster un- der conditions gi ven in Theorem 4, it is natural to in vesti- gate which factors contribute to the conv ergence speed, i.e., λ D +2 ( Ω ) , giv en that each cluster is locally optimized for large λ D +2 ( Ω c ) . Theorem 5 ( Ensemble con verg ence speed λ D +2 ( Ω ) ) . If the conditions in Theorem 4 are met for collecti ve motions, the con vergence-indicating eigen value of the ensemble stress matrix, i.e., λ D +2 ( Ω ) for Ω in (30), is upper bounded by λ D +2 ( Ω ) ≤ min c ∈{ 1 ,...,C } ( λ D +2 ( Ω c ) + ρ c ) . (35) Here, defining v c ∈ R N c as the eigen vector associated with λ D +2 ( Ω c ) , ρ c = β P C k  = c   v c |B ck   2 2 with β in the spectral norm constraint in P 2 , and v c |B ck the subv ector of v c corre- sponding to the bridging nodes to cluster k . Proof . The smallest non-zero eigenv alue can be defined as λ D +2 ( Ω ) = min v ⊥ Null( Ω ) ∥ v ∥ =1 v ⊤ Ω v ≤ v ⊤ Ω v v ⊤ v (36) for any candidate vector v ∈ R N , with v ⊥ Null( Ω ) . Consider the candidate vector ˜ v c =  0 ⊤ v ⊤ c 0 ⊤  ⊤ ∈ R N , which is the eigen vector v c ∈ R N c associated with λ D +2 ( Ω c ) , correctly zero-padded to length N . Lemma 2 in Appendix A shows that ˜ v c ⊥ Null( Ω ) and ∥ ˜ v c ∥ = 1 , and thus ˜ v c satisfies the bound (36), i.e., λ D +2 ( Ω ) ≤ ˜ v ⊤ c Ω ˜ v c . Further extending this bound, we obtain λ D +2 ( Ω ) ≤ ˜ v ⊤ c Ω ˜ v c = ˜ v ⊤ c   ˜ Ω c + C X k  = c ˜ Ω k   ˜ v c (37a) = ˜ v ⊤ c ˜ Ω c ˜ v c + ˜ v ⊤ c ( C X k  = c ˜ Ω k ) ˜ v c , (37b) with ˜ v ⊤ c ˜ Ω c ˜ v c = v ⊤ c Ω c v c = λ D +2 ( Ω c ) (38a) ˜ v ⊤ c ( C X k  = c ˜ Ω k ) ˜ v c = C X k  = c v ⊤ c |B ck Ω k |B ck v c |B ck , (38b) Bridging Nodes: (a) Bridging nodes span 2D (b) Rank-deficient bridge (c) Rank-deficient bridge Fig. 5. An illustration of local flexibility introduced by rank-deficient bridging nodes. The shaded area outlines the desired configuration and the dashed geometry sho ws the potential ambiguity . where v c |B ck denotes the sub-vector of v c corresponding to the bridging nodes B ck shared between cluster c and k , and Ω k |B ck is the subblock of Ω k corresponding to those bridging nodes. W e could further relax (38b) using v ⊤ c |B ck Ω k |B ck v c |B ck ≤ λ max  Ω k |B ck    v c |B ck   2 2 (39a) ≤ λ max ( Ω k )   v c |B ck   2 2 (39b) ≤ β   v c |B ck   2 2 , (39c) where (39b) is the relaxation by the Cauchy interlacing theo- rem [31] which upper bounds the maximum eigenv alue of a principal submatrix by that of the parent matrix, and (39c) is spectral norm constraint of P 2 in the stress design. Substituting (39c) back to (38b) giv es us ρ c = β P C k  = c   v c |B ck   2 2 . Finally , since this inequality (37) is satisfied for all c = 1 , ..., C , we adopt the tightest one, i.e., bound (35). ■ W ith the upper bound (35) we can infer that the con vergence speed of the ensemble graph is bottlenecked by two factors: the slowest cluster indicated by λ D +2 ( Ω c ) , and the coupling strength of the bridging nodes indicated by ρ c . This provides some guidelines for the design of the clusters, i.e., to make sure each cluster conv erges fast using the proposed algorithm P 2 , and to ensure the bridging nodes connect the clusters well. The design of the optimal bridging remains an open problem; howe ver , heuristics such as using more bridging nodes for a fixed β are proven to be effecti ve in practice. D. Extension to P artial-collective Behaviors Follo wing the conclusion of the pre vious subsection, if the bridging nodes are not selected such that they span D - dimensional space, the clusters may not be aligned in motion. Instead, the swarm can be intentionally designed to exhibit partial-collectiv e motions, allowing for cluster-wise flexibility . In the most extreme case, if no bridging nodes are present, the clusters act as completely independent systems. T o maneuver an affine formation, the global ensemble fundamentally requires a base set of at least D + 1 affinely independent leaders to span the space [17]. Howe ver , when the swarm features rank-deficient bridges, this global leader set alone is insufficient to uniquely determine the configuration of e very sub-cluster . The bridging nodes V B c within cluster c act as virtual leaders, providing constraints dictated by the neighboring cluster . T o guarantee that cluster c is affinely localizable, meaning its local shape and affine transformations 10 T ABLE I B E N C H M A R K R E S U LT S O N R U N T I M E ( I N [ S E C O N D S ] ) test case Random-5 Random-8 Random-50 Random-100 Random-500 Circular-10 Cuboctahedron T runc. Icosa. Letter-“W” N 5 8 50 100 500 10 12 60 260 D 2 2 2 2 2 2 3 3 3 proposed ( P 2 ) 6.7e-3 10.3e-3 4.8 N/A N/A 13.8e-3 16.3e-3 7.7 N/A proposed-usi ( P 3 ) N/A N/A N/A N/A N/A 8.9e-3 8.6e-3 0.27 N/A proposed-m.c. 1 – – – 17.6 88.5 – – – 18.5 2 Lin et al. 4.1e-3 7.1e-3 N/A N/A N/A 9.5e-3 N/A N/A N/A Y ang et al. 0.080e-3 0.10e-3 0.55e-3 0.97e-3 4.8e-3 0.2e-3 0.1e-3 0.6e-3 2.5e-3 Xiao et al. 519.0e-3 108.3 N/A N/A N/A 322.7 N/A N/A N/A Li et al. 0.307e-3 0.56e-3 2.7e-3 5.3e-3 31.1e-3 0.65e-3 0.7e-3 3.7e-3 15.2e-3 “N/A”: Not applicable due to memory constraints, impracticably long runtime, or inv alid results. “–”: Multicluster operations not needed. 1 The proposed-m.c. is configured with 2 clusters of 60 nodes for “Random-100”, 10 clusters of 60 nodes for “Random-500”, and 4 clusters of 58 nodes for “Letter-W”. 2 The computation is only needed for 1 cluster, as all clusters are identical up to affine transformations. (a) Circular-10 in 2D (b) Cuboctahedron in 3D (c) T runcated-Icosahedron in 3D (d) Letter-“W” in 3D Fig. 6. The special configurations used for the benchmark. For 3D configu- rations, the edges only outline the geometry for clear visualization. are fully determined and controllable, the actual leaders within the cluster must compensate for the bridge’ s rank deficiency . They must jointly satisfy the local condition rank  ¯ P (( V l ∩ V c ) ∪ V B c )  = D + 1 , ∀ c = 1 , ..., C (40) where V l is the global subset of chosen leaders, and V B c is the set of bridging nodes attached to cluster c . Physically , (40) dictates how to control the unconstrained degrees of freedom with leaders introduced by a rank-deficient bridge. Se veral 2D examples are visualized in Fig. 5, where two clusters with different bridging node setups are considered. In Fig. 5 (a), the bridging nodes satisfy conditions for collectiv e behaviors, thus no additional leaders are required beyond the 3 for global maneuver . In (b), the bridging nodes form a line, introducing loose local degrees of freedom, as can be seen by the dashed outline. In this case, at least 1 leader should be allocated to each cluster . Similarly , (c) has 1 bridging node and requires 2 leaders for each cluster to ensure global localizability . V I . S I M U L A T I O N S In this section, we numerically validate the proposed al- gorithms and the multicluster control frame work. W e first perform a benchmark comparing with multiple state-of-the- art stress matrix design algorithms on sev eral performance metrics. W e then demonstrate the conv ergence and the inter- cluster flexibility of our proposed multicluster control. All code and plots are publicly av ailable 1 . A. Benchmarks of Stress Design The algorithms are benchmarked on a Linux platform with an AMD Ryzen 9 7900X processor and 64GB of DRAM. W e select some representative state-of-the-art candidates from all existing methods, including Lin et al. [16], Xiao et al. [26], Y ang et al. [23], and Li et al. [24]. The proposed solutions include the generic stress design method P 2 from Section III with 3 sets of hyperparameters (namely γ = 0 . 1 , β = 1 , and α ∈ { 0 . 5 , 1 . 5 , 5 } ), the unique stress identifier (USI) method P 3 from Section IV for reduced complexity , and the multicluster division from Section V. The optimization-based algorithms are solved using MOSEK [32] with the CVXPY interface [33], and all algorithms perform 100 Monte Carlo runs for each case. W e test our solutions using random con- figurations of different sizes, namely N ∈ { 5 , 8 , 50 , 100 , 500 } with a naming conv ention such as “Random-5”, “Random-8”, etc. Some rotationally symmetric configurations are considered as well as shown in Fig. 6, including “Circular-10”, which is a 10 -node regular decagon, a 12 -node “Cuboctahedron”, and a 60 -node “Truncated Icosahedron”. Lastly , we use a customized multicluster configuration “Letter-W” in Fig. 6(d) to demonstrate the partial-collecti ve behaviors of the multi- cluster control, where the 4 segments of the configurations are identical up to affine transformations. Our performance metrics include (a) the average runtime of the algorithms, which indicates their complexity and scalabil- ity; (b) the av erage degree 2 M / N of the resulting graph, which reflects the sparsity and the general communication ov erhead; (c) eigen value λ D +2 ( Ω ) , which indicates the conv ergence speed of the formation; and (d) the spectral ef ficiency η = λ D +2 ( Ω ) N 2 λ N ( Ω ) M , (41) which assesses the con vergence speed, given the same commu- nication capacity . For the multicluster setup, we ev aluate the equiv alent stress matrix Ω = P C c =1 ˜ Ω c in the mean dynamics. The results of the runtime are reported in T able I and statistics of other metrics are shown in Fig. 7. 1 https://github .com/asil- lab/zli- large- afc 11 W e first evaluate the complexity of the algorithms through the a verage runtime in T able I, from which we generally observe that the optimization-based algorithms consume a longer time to compute, compared with the analytical Y ang et al. and Li et al. In particular , Lin et al. is constrained by the (manual) prescription of the graph, and Xiao et al. (MISDP) is only capable of small-scale networks due to the NP-hard complexity . The generic version of the proposed algorithm is capable of handling medium-sized networks ( N < 100 ), abov e which the multicluster division results in a reasonable runtime. In the case of rotationally symmetric configurations, such as circular and cuboctahedron, the USI method with P 3 effecti vely reduces the runtime compared to the generic solu- tion P 2 . In case of the truncated icosahedron, the USI method giv es a roughly 30 × speedup. For special geometries such as Fig. 6(d), where a subconfiguration is repeated multiple times up to affine transformations, it is natural to divide them into sev eral clusters, for which the stress matrices also repeat and only need to be computed once. The sparsity of the graph is shown in the top figure in Fig. 7. Existing state-of-the-art methods typically yield sparse solutions, whereas the proposed approach allows explicit con- trol ov er sparsity through the hyperparameter α . Particularly , for small values (i.e., α = 0 . 5 , which is around the critical value α ∗ introduced in Section III-C), it produces graphs that are comparably sparse to those obtained by state-of-the-art methods. F or a large α = 5 , it approaches the upper bound cor- responding to a complete graph. Note that in the multicluster implementations of “Random-100” and “Random-500”, this upper bound cannot be reached ev en for large α , since non- bridging nodes among the clusters do not interact. This feature naturally promotes a sparse solution and can further be adapted to accommodate, e.g., limits of communication ranges. The conv ergence speed and spectral efficienc y are also ev aluated in Fig. 7. From the smallest non-zero eigen value λ D +2 ( Ω ) , we observe that the analytical algorithms of Y ang et al. and Li et al. yield v alues that are 1 or 2 orders of magnitude smaller for medium- to large-scale networks. This may be considered a trade-off of their sparse nature. The optimization-based algorithms, including the proposed, generally yield larger λ D +2 ( Ω ) , meaning a faster con vergence speed since systematic searches for optimality are in volv ed. Note that for the proposed algorithm with α = 0 . 5 , the con ver gence is much faster than the analytical solutions while achieving similar levels of sparsity across the various cases. This is also seen from the spectral ef ficiency plot in the bottom figure of Fig. 7, where the proposed algorithms better exploit the communication capacity into accelerating the conv ergence more efficiently . W ith α tuned larger , the resulting graph becomes denser while further accelerating the con vergence. Howe ver , larger α may be less fav ored in practice as the con- ver gence boost is increasingly marginal, while compromising the graph sparsity . T o summarize the ev aluation on stress design, our proposed algorithms have sparse results, fast con vergence, and high spectral ef ficiency while admitting acceptable complexity for medium and large-scale networks. Additionally , the y are robust to generic and nongeneric configurations in 2D and 3D, and Average Degree (2 M/ N ) λ D +2 ( Ω ) Configuratio ns (Nodes) Spectral E ffi ciency ( η ) proposed proposed-m.c. α =0 . 5 α =1 . 5 proposed proposed-m.c. α =5 proposed proposed-m.c. Random (5) Random (8) Circular (10) Cuboctahed ron (12) Random (50) Trunc . Icosahedron (60) Random (100) Random (500) Y ang et al. Xiao et al. Lin et al. Li et al. upper b ound Fig. 7. Benchmark performances of the proposed algorithms and the state- of-the-arts. The plotted metrics include the graph sparsity characterized by av erage degree of the graph 2 M / N (upper), the conv ergence-speed-indicating eigen value λ D +2 ( Ω ) (middle), and the spectral efficienc y (41) (lower). yield consistently superior results as compared to state of the art across various benchmarks. B. Multicluster Contr ol with Collective Motions W e no w validate the proposed multicluster control frame- work using the “Random-100” test case. Recall from Section V that a certain number of bridging nodes are required for the clusters to engage in collectiv e motions, which also affects the conv ergence speed in the mean dynamics. Fig. 8(a) demonstrates an example of 2 clusters with dif ferent numbers of bridging nodes, and (b) shows the con ver gence speed through the eigen value λ D +2 ( Ω ) of the cluster ensemble in the left plot, and the tracking error conv ergence on the right. Since the “Random-100” is a 2D configuration, the mini- mum number of bridging nodes needed is D + 1 = 3 and they must span the 2D space. It can be clearly seen from Fig. 8(b) that λ D +2 ( Ω ) is 0 and the control tracking error does not decrease when there are only 2 bridging nodes. This shows the lack of collecti ve behaviors due to the loose degrees of freedom caused by insufficient bridging nodes. As we increase the bridging nodes, we observe larger λ D +2 ( Ω ) , which suggests a faster con ver gence. Correspondingly , the tracking error presents steeper conv ergence curves with more bridging nodes. 12 Bridging No des: (a) ”Random-100” configuration with 2 clusters. The number of bridging nodes is varied among 2 (not sufficient for collective motion), 5, 10, and 20. Index Time [s] λ D +2 ( Ω ) 20 10 5 2 Number of Bridging Nodes (b) The conver gence speed through the eigen values (left, only plotted the first 10) and the control tracking error (right). Fig. 8. The results of the multicluster control with collective affine motions. The number of bridging nodes is varied to validate the claim about the upper bound on the ensemble graph’ s eigen values. C. Swarm with Inter-Cluster Flexibility W e finally demonstrate the feature of inter-cluster flexibility of the multicluster control framew ork using the “Letter-W” configuration as shown in Fig. 6(d). The configuration is constructed so that each letter segment (considered a cluster) is identical up to an affine transformation, allowing the stress matrix to be computed for only a single segment. T o allow intercluster flexibility , the clusters are interfaced with 4 nodes on a plane, and the remaining de grees of freedom are tightened by leaders within each cluster . Then by applying local affine transformations, the clusters may present other geometries than the “Letter-W” by simply maneuvering the leader positions. Fig. 9 demonstrates the maneuvering patterns, where all nodes are initialized randomly in 3D, and they sequentially con verge to the target formation of a “Bar”” shape, followed by a “Letter-V” shape, and finally the “Letter-W” while admitting a collectiv e translation movement. This validates the feature of inter -cluster flexibility with partial collecti ve motions of the proposed multicluster control framework. V I I . C O N C L U S I O N In this work, we propose a con vex formulation for de- signing the stress matrix used for affine formation control, aiming to optimize the conv ergence speed while keeping the network sparse. For the rotationally symmetric configurations, we show that the optimal solution is highly structured, in which the stress values are identical for equiv alent edges, and thus propose a reduced formulation for lower computa- tion complexity . Compared with the state-of-the-art solutions through a comprehensiv e benchmark, we conclude that our proposed formulation is superior in the con vergence speed and sparsity while maintaining acceptable computation complexity . Having the optimal design of the network, we also propose a multicluster control framework in which large-scale networks are di vided into smaller networks with the ov erlapping nodes t = 0s t = 20s t = 40s t = 60s Fig. 9. Large-scale formation maneuvers with inter-cluster flexibility . The swarm is initialized at a random configuration in 3D. Multicluster control is deployed to achiev e the “Letter-W” and sev eral other configurations using the same optimized stress matrix. configured to present both collectiv e and noncollective beha v- iors. In future work, to realize large swarm implementations, we aim to consider large-scale collision av oidance, explicit modeling and optimization of the communication cost. As a direct application of the multicluster framew ork, extensions of the stress-based control used in this work can be used for relev ant scenarios. Furthermore, the selection of the bridging nodes can also be optimized for collective motions. A P P E N D I X A P R O O F S O F P R O P E RT I E S A N D L E M M A S Property 1 ( P ermutation in variance of tr( Ω ) ) . The trace of the stress matrix is in variant under permutations, i.e., tr ( Ω ) = tr  Π k ΩΠ ⊤ k  Proof . Using the cyclic property of the trace operator , tr  Π k ΩΠ ⊤ k  = tr  Π ⊤ k Π k Ω  = tr( Ω ) since permutation matrices are orthogonal, i.e., Π ⊤ k Π k = I N . ■ Proof of Lemma 1 . Proof . W e first show the inv ariance property of the objective function. It is straightforward to see that   ¯ Π E ,k ¯ ω   1 = ∥ ¯ ω ∥ 1 for any permutation ¯ Π E ,k ∈ { 0 , 1 } ¯ M × ¯ M since permuting the elements of ¯ ω does not change the L1 norm. The trace of the stress matrix tr( Ω ) is also permutation in v ariant, sho wn in Property 1 in Appendix A. W e could thus conclude that the objectiv e (7a) satisfies permutation in variance in Definition 3. Next, we show the in variance property of the feasible set. Observe from (13) that the permutations of the stress matrix Π k ΩΠ k ⊤ is a similarity transforms of Ω , which preserves the eigen values. As such, Constraints (7c) and (7d) that limit the range of the eigen values still hold under permutations. As for the equality constraint (7e), we extend the definition (12) to ¯ R k ¯ P = ¯ P Π k where ¯ R k = diag ( R k , 1) ∈ R ( D +1) × ( D +1) . Observe that ¯ B diag  ¯ Π E ,k ¯ ω  ¯ B ⊤ ¯ P ⊤ (13) = Π k ΩΠ k ⊤ ¯ P ⊤ (42a) = Π k Ω ¯ P ⊤ ¯ R ⊤ k (42b) (7e) = 0 , (42c) which prov es the permutation in variance of constraint (7e). As such, problem P 1 is permutation-in variant. ■ 13 Lemma 2 ( Orthogonality of ˜ v c and Null( Ω ) ) . Consider ˜ v c =  0 ⊤ v ⊤ c 0 ⊤  ⊤ ∈ R N , which contains the eigen vector v c ∈ R N c associated with λ D +2 ( Ω c ) , zero-padded to the length of N . Then ˜ v c ⊥ Null( Ω ) if rank( Ω ) = N − D − 1 . Proof . If rank( Ω ) = N − D − 1 , then augmented configuration matrix ¯ P ( V ) fully span Null( Ω ) . Then we could verify that ˜ v ⊤ c 1 N = v ⊤ c 1 N c = 0 since vector 1 N c is in the nullspace of Ω c and is thus orthogonal to v c . Similarly , we could verify that ˜ v ⊤ c P ⊤ = v ⊤ c P ⊤ c = 0 ⊤ D . Therefore, we could establish that ˜ v c ⊥ Null( Ω ) and ˜ v ⊤ c ˜ v c = 1 since ˜ v c is zero-padded. ■ R E F E R E N C E S [1] Z. Li, G. Leus, and R. T . Rajan, “Fast multiagent formation stabilization with sparse uni versally rigid frame works, ” in 2025 33rd Eur opean Signal Pr ocessing Confer ence (EUSIPCO) . IEEE, 2025, pp. 2392–2396. [2] W . Zhu, S. O ˘ guz, M. K. Heinrich, M. Allwright, M. W ahby , A. L. Chris- tensen, E. Garone, and M. Dorigo, “Self-organizing nervous systems for robot swarms, ” Science Robotics , vol. 9, no. 96, p. eadl5161, 2024. [3] G. Sun, R. Zhou, Z. Ma, Y . Li, R. Groß, Z. Chen, and S. Zhao, “Mean-shift exploration in shape assembly of robot swarms, ” Natur e Communications , vol. 14, no. 1, p. 3476, 2023. [4] S. Na, T . Rou ˇ cek, J. Ulrich, J. Pikman, T . Krajn ´ ık, B. Lennox, and F . Arvin, “Federated reinforcement learning for collective navigation of robotic swarms, ” IEEE T ransactions on cognitive and developmental systems , vol. 15, no. 4, pp. 2122–2131, 2023. [5] L. Bartolomei, L. T eixeira, and M. Chli, “Fast multi-uav decentralized exploration of forests, ” IEEE Robotics and Automation Letters , vol. 8, no. 9, pp. 5576–5583, 2023. [6] M. Jurt, E. Milner, M. Sooriyabandara, and S. Hauert, “Collectiv e transport of arbitrarily shaped objects using robot swarms, ” Artificial Life and Robotics , vol. 27, no. 2, pp. 365–372, 2022. [7] G. Zhao, H. Cui, and C. Hua, “Hybrid event-triggered bipartite consen- sus control of multiagent systems and application to satellite formation, ” IEEE T ransactions on Automation Science and Engineering , vol. 20, no. 3, pp. 1760–1771, 2023. [8] K.-K. Oh, M.-C. Park, and H.-S. Ahn, “ A survey of multi-agent formation control, ” Automatica , vol. 53, pp. 424–440, 2015. [9] Z. Li and R. T . Rajan, “Geometry-aw are edge-state tracking for robust affine formation control, ” IEEE Open Journal of Contr ol Systems , 2026. [10] G.-P . Liu and S. Zhang, “ A survey on formation control of small satellites, ” Pr oceedings of the IEEE , vol. 106, no. 3, pp. 440–457, 2018. [11] R. Babazadeh and R. R. Selmic, “Directed distance-based formation con- trol of nonlinear heterogeneous agents in 3-d space, ” IEEE T ransactions on Aerospace and Electronic Systems , vol. 59, no. 3, pp. 3405–3415, 2022. [12] X. Fang and L. Xie, “Distributed formation maneuver control using complex laplacian, ” IEEE T ransactions on Automatic Control , vol. 69, no. 3, pp. 1850–1857, 2023. [13] H. Cheng and J. Huang, “ A general framework for the bearing-based formation control, ” IEEE T ransactions on Automatic Contr ol , vol. 70, no. 6, pp. 3603–3616, 2024. [14] J. Zheng and A. Jamalipour, W ir eless sensor networks: a networking perspective . John W iley & Sons, 2009. [15] R. Heusdens and G. Zhang, “Distributed optimisation with linear equal- ity and inequality constraints using pdmm, ” IEEE T ransactions on Signal and Information Processing over Networks , 2024. [16] Z. Lin, L. W ang, Z. Chen, M. Fu, and Z. Han, “Necessary and sufficient graphical conditions for affine formation control, ” IEEE Tr ansactions on Automatic Control , vol. 61, no. 10, pp. 2877–2891, 2015. [17] S. Zhao, “ Affine formation maneuver control of multiagent systems, ” IEEE T ransactions on Automatic Control , vol. 63, no. 12, pp. 4140– 4155, 2018. [18] Q. Y ang, X. Zhang, H. Fang, M. Cao, and J. Chen, “Joint estimation and planar affine formation control with displacement measurements, ” IEEE Tr ansactions on Control Systems T echnology , vol. 33, no. 1, pp. 92–105, 2024. [19] R. Connelly and S. D. Guest, Fr ameworks, tense grities, and symmetry . Cambridge Uni versity Press, 2022. [20] Y . Kim and M. Mesbahi, “On maximizing the second smallest eigen value of a state-dependent graph laplacian, ” IEEE T ransactions on Automatic Contr ol , vol. 51, no. 1, pp. 116–120, 2006. [21] R. Olfati-Saber, J. A. Fax, and R. M. Murray , “Consensus and coop- eration in networked multi-agent systems, ” Pr oceedings of the IEEE , vol. 95, no. 1, pp. 215–233, 2007. [22] S. D. Kelly and A. Micheletti, “ A class of minimal generically univer- sally rigid frameworks, ” arXiv pr eprint arXiv:1412.3436 , 2014. [23] Q. Y ang, M. Cao, H. Fang, and J. Chen, “Constructing univ ersally rigid tensegrity framew orks with application in multiagent formation control, ” IEEE T ransactions on Automatic Control , v ol. 64, no. 1, pp. 381–388, 2018. [24] H. Li, H. Chen, X. W ang, Z. Li, and L. Shen, “Distributed framework construction for affine formation control, ” IEEE T ransactions on Auto- matic Contr ol , 2025. [25] Q. Y ang, M. Cao, and B. D. Anderson, “Growing super stable tensegrity framew orks, ” IEEE T ransactions on Cybernetics , vol. 49, no. 7, pp. 2524–2535, 2018. [26] F . Xiao, Q. Y ang, X. Zhao, and H. Fang, “ A framework for optimized topology design and leader selection in affine formation control, ” IEEE Robotics and Automation Letters , vol. 7, no. 4, pp. 8627–8634, 2022. [27] R. Connelly , “Generic global rigidity , ” Discr ete & Computational Ge- ometry , vol. 33, pp. 549–563, 2005. [28] S. J. Gortler and D. P . Thurston, “Characterizing the uni versal rigidity of generic frameworks, ” Discrete & Computational Geometry , vol. 51, no. 4, pp. 1017–1036, 2014. [29] Y . Xu, S. Zhao, D. Luo, and Y . Y ou, “ Affine formation maneuver control of linear multi-agent systems with undirected interaction graphs, ” in 2018 IEEE confer ence on decision and control (CDC) . IEEE, 2018, pp. 502–507. [30] E. De Klerk, Aspects of semidefinite pro gramming: interior point algo- rithms and selected applications . Springer, 2002. [31] S.-G. Hwang, “Cauchy’ s interlace theorem for eigenv alues of hermitian matrices, ” The American mathematical monthly , vol. 111, no. 2, pp. 157–159, 2004. [32] M. ApS, “Mosek modeling cookbook, ” 2020. [33] S. Diamond and S. Boyd, “CVXPY: A Python-embedded modeling lan- guage for conve x optimization, ” Journal of Machine Learning Researc h , vol. 17, no. 83, pp. 1–5, 2016.