Hierarchic distributed stabilization of a class of three-dimensional formations for underactuated agents

The paper presents a distributed hierarchic control design solving a formations problem in three dimension. The agents are modelled as thrust-underactuated rigid bodies. The class of formations addressed encompasses path following of closed convex paths with shape, size, orientation, relative displacements and relative agents’ on-path angular positions all seen as formation parameters. This problem can be seen as a generalization of circular formation problems which acquired particular attention in the ﬁeld of multi-agents, and the agents model used is one that usually models certain unmanned aerial vehicles and other autonomous vehicles. The solution relies on reduction-based set stabilization where the problem is broken down into simpler nested sub-problems. The results illustrate how addressing complex control problems using hierarchical set stabilization is natural, simpli-ﬁes the solution, and provides useful properties. The results are illustrated via simulations.


INTRODUCTION
Formation control continues to receive wide attention in control theory and practice. Recent advances in computation, communications and sensor technologies provide the capability to build teams of relatively small and inexpensive robots. These can perform cooperative tasks in air, on land, underwater or in space; and often rely on formation control to fulfill these tasks. The span of applications for such systems is quite diverse [1]. Examples include: environmental monitoring and exploration, surveillance, search and rescue, reconnaissance, space interferometry, and mobile and reconfigurable sensors [1].
Numerous formation control methods have been introduced in literature [2]. This paper utilizes an approach that has not received particular interest [3,4]. The approach is hierarchical in nature and relies on theories of set stabilization and reduction [5]. Simply put, the desired formation behaviour is first expressed as a goal set in the agents state space. A hierarchy then ensues by breaking down the problem into l ∈ ℤ + sub-problems of stabilizing nested sets Γ 1 ⊃ ⋯ ⊃ Γ l , with Γ l encoding the formation (goal set), and sub-problem i is that of stabilizing Γ i relative to Γ i−1 (Γ 0 being the state space). The motivation and contribution aspects of the paper can be seen as threefold.
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. First, the paper illustrates the potential benefit of hierarchical set stabilization in formation control problems. Formations are usually given by certain relations between the agents; in addition the agents are required to attain certain dynamics while in formation. Being the generalization of equilibrium stabilization set stabilization presents a natural framework to consolidate several specifications into one: the goal set. Set stabilization also naturally prompts a hierarchic approach where subsets of the specifications provide nested sets containing the goal set. Hierarchy is favorable in solving control problems from a number of perspectives. From a design perspective, with the appropriate tools, a complex problem can be solved more easily by breaking it down into simpler decoupled sub-problems. Individual subproblems can be of independent interest. From a practical perspective, if, with initial design, the system satisfies specifications that are higher in the hierarchy, then a hierarchic design would stabilize the system to lower specifications without violating the higher ones.
Hierarchy has been used in different contexts in formation control, see, for instance, [6]. See also reduction-based results such as [7,8]. Combining hierarchy and set stabilization, as done here, helps to elucidate the stability analysis, especially for goal behaviours expressed by unbounded sets, a common instance in formation control.
Second, the formation problem addressed here can be seen as a generalization of circular formations problems which received significant attention in literature. The breadth of applications is quite wide [9]. Examples include patrolling and surveillance [10,11], sensing and data collection [12], target tracking and source seeking [13], search and rescue [14], environmental monitoring [15], and aerospace applications [16,17]. The generalization presented here provides certain advantages which makes the problem more oriented to the practical setting. Refer to Section 2.1 for a discussion on the problem potential. In addition, the hierarchic solution approach provides extra advantages on its own. Refer to the discussion at the end of Section 3.
Numerous methods and different instances of the circular formation problem can be found in literature. Common approaches include optimization methods [18], virtual structure control [], model predictive control [20], artificial potential functions [21], leader follower approach [22], and game theory [23]. Cyclic pursuit was studied in [24], Lyapunov guidance vector field was used in [25] for circular orbit stabilization for unmanned aerial vehicles (UAVs), and in [26] circular formations were stabilized using modified Kuramoto model. The problem addressed here is generalized to include several instances: co-path and different paths, co-formation plane and different planes, and uniform and non-uniform formations. The approach further breaks the problem down into sub-problems of potential independent interest. Refer to Sections 2.1 and 3.
Regarding the control algorithm used, the main difference between it and previous methods lies more in the conceptional design (hierarchic set stabilizing) than in the particulars of the controllers. The main focus here is on the breakdown of the problem to simpler sub-problems that can be solved relatively simpler than by using a direct approach, and that could allow for the use different methods to solve the particular sub-problems.
The approach of the paper has been used to solve for formation control in [3] and [4]. Those results can be seen as precursors to this work. In [3], a circular formation problem was solved for dynamic unicycles using hierarchical set stabilization. A group of unicycles were required to follow a common circular path with prescribed radius, and desired on-path separations. The problem was solved using three hierarchical sub-problems. In [4], a circular formation problem was solved for steered kinematic particles in three-dimensional space. The particles did not necessarily need to follow a common path. The formation was determined by the spacings between the different paths, and the relative angular positions of the particles on their paths. These aspects are shared with the problem addressed here. The solution in [4] broke down the problem into three sub-problems where the last is similar to that of [3]. Although the problem addressed here shares a number of aspects with that of [4], yet it is considered more amenable to practical implementation where the agents are dynamic and underactuated, the paths are not necessarily circular and the communication graph is directed (in [4] it was undirected). The results here, as well as in [3] and [4], utilize the globally stabilizing feedbacks for path following of strictly convex paths developed in [27], and the reduction theorem developed in [5] is used to consolidate the sub-problems solutions to achieve asymptotic stabilization of the goal set.
The third contribution aspect of the paper concerns the agents model used. Different models have been considered in literature when studying multi-agent formations. These include mobile robots and UAVs [6,10,11,13,18], spacecrafts and satellites [21,28,29], and surface and underwater vehicles [30]. The agents model used in this paper can be used for underwater, and aerial or space flying robots. However, that model would usually be accompanied with force fields, and/or with coupled degrees of actuation. While such forces and couplings are not considered here it is expected that the hierarchy could allow to address such challenges. Refer to Section 7 for a discussion on how to extend the current results to address more practical agent models with the previously mentioned limitations.
A particular point of interest of the model used here is its underactuations. Handling actuation limitations is a significant challenge for multi-agent systems and problems. The paper shows how the approach presented can be beneficial in addressing this challenge. Instead of developing a particular solution for the formation problem in question, underactuation is addressed by controlling the attitudes to an orientation where the unactuated degrees are not needed. This is possible by achieving a decoupling (stemming from the hierarchy) between the underactuation aspect of the problem and the formation ones. The result is a relatively simple solution that can be utilized in other problems, and that can incorporate available algorithms for attitude control. Underactuated rigid bodies in ℝ 3 has been used as model for many autonomous flying vehicles. For instance, quadrotor helicopter and vertical take-off and landing (VTOL) aircraft can be modelled as degree-two underactuated rigid bodies with one thrust and three torques. Control for such autonomous systems has received considerable attention [31][32][33][34]. The design in this paper utilizes the general idea of using attitude stabilization to acquire desired thrust orientation. However, the problem here is more involved with different specifications, and information exchange constraints.
Notation. For scalars a 1 , … , a n , diag(a 1 , … , a n ) denotes the diagonal matrix with entries a i . If A, B are two matrices, diag(A, B) denotes the block diagonal matrix with blocks A, B, and A ⊗ B denotes their Kronecker product. For vectors a, b, col(a, b) denotes the concatenation of a and b vertically. The index set {1, … , n} is denoted by n, and the n-vector of ones is denoted by 1. For x ∈ ℝ, x mod 2 denotes its value modulo 2 , and if , x ∈ ℝ then x = mod 2 states that (ℝ) will be used to denote the class of C 1 saturation functions : ℝ → ℝ such that for all y ∈ ℝ, (y)y ≥ 0, | (y)y| < 1, and (y) = 0 only when y = 0. Without loss of generality, it will be assumed here that ′ (0) = 1. For x ∈ ℝ 3 , x × denotes its skew symmetric matrix representation, where x × y = x × y for any y ∈ ℝ 3 . For a dynamical system Σ :ẋ = f (x), (t, x 0 ) denotes the solution with initial condition x 0 at t = 0. Finally, for a closed set Γ, ‖ ⋅ ‖ Γ denotes the point to set distance, and for > 0, B (Γ) := {y ∈  : ‖y‖ Γ < }.

PROBLEM FORMULATION
A system of n ≥ 2 rigid bodies in three-dimensional Euclidean space is considered. Figure 1 depicts agent i.
denote the body frame in inertial frame  = span( 1 , 2 , 3 ) (the standard basis). It is assumed here that each rigid body is provided with two thrusts f i The states of agent i are given as: x i ,ẋ i ∈ ℝ 3 the centre of mass and linear velocity in , the attitude matrix, and w i ∈ ℝ 3 the angular velocity in  i . The state vector of agent i is given by 3 , and the collective state vector is given by = ( 1 , … , n ) ∈  :=  1 × ⋯ ×  n . The system model is given bÿ where ⊤ is input thrust, and m i is the mass of rigid body i. Let x = (x 1 , … , x n ),ẋ = (ẋ 1 , … ,ẋ n ), B = (B 1 , … , B n ) and b 1 = (b 1 1 , … , b n 1 ).

Assumption 1.
The angular dynamics of (1) are neglected, that is, it is assumed that desired angular velocities w i can be achieved instantaneously.
Remark 1. Due to hierarchic nature of the approach, Assumption 1 can be relaxed to: desired angular velocities can be achieved exponentially, regardless of the form of torques used or their degree of actuation. By this, and on account of available literature on angular dynamics stabilization, the instantaneous assumption is used.

Assumption 2.
The information exchange between the rigid bodies is modelled by a directed sensor graph  where • agent i has access to its orientation B i ; its linear and angular velocitiesẋ i , w i ; and • relative positions, orientations, linear and angular velocities expressed in  i of agents visible to it per .  will be considered to be static and balanced, with a globally reachable node [35]. N (i ) will be used to denote the set of nodes connected to i according to . L will denote the Laplacian of  with L i its ith row, and L (3) The following lemma will be used.

Lemma 1.
If L is the Laplacian of an n-node digraph with a globally reachable node, then for sufficiently small k > 0, has eigenvalues with negative real parts except for a simple eigenvalue at 0.
The proof is provided in Appendix A.1.
Part of the formation problem addressed here entails path following of a closed convex path. Consider a smooth Jordan curve ℭ ∈ ℝ 2 that is strictly convex. From lemma 3.1 in [27], there exists a regular parametrization = ( 1 , 2 ) : S 1 → ℝ 2 of ℭ such that, for each ∈ S 1 , the angle of the tangent vector ′ ( ) is mod 2 , refer to Figure 2(a). Using a frame  with axes p 1 , p 2 , p 3 , a strictly convex path in ℝ 3 can be given by where P = [p 1 p 2 p 3 ], and y c ∈ ℝ 3 denotes the path centre. The frame  determines the orientation of the path in ℝ 3 where p 1 , p 2 span the path plane, refer to Figure 2(b). This path will be denoted ( ,  ). Now, suppose that agent i is following the path with forward heading b i 1 , and let ) . ( By this, i 1 equals the angle of the tangent to the path at x i , and i determines the on-path angular position of rigid body i. Let

Problem statement
For system (1) of n rigid bodies with digraph  modelling the information exchange between them. The Convex-Path-Following Formations problem (C x PFF) is defined as the distributed stabilization problem of the goal set where ∈ ℝ 3n , ∈ ℝ n ,ū :  i → ℝ + a C 1 function, and This problem is a combination of the following specifications.
(i) Convex path following: in C x PFF each agent follows, the same or different, strictly convex path ( ,  ) with stationary centre (dependent on initial conditions). The path plane orientation is determined by p 3 , that is, if the agent heading is in this plane then p 3 ⋅ b i 1 = 0 1 . Also, while path following, the heading should be in the direction of the linear velocity, that is,ẋ i ⋅ b i 1 equals the forward speed. Let this beū( i ) > 0. By proposition 3.4 in [27], in order for the agent to stay on the path, the angular velocity on the path, that is, the component of w i orthogonal to the path should bēu . By this, specification (i) can thus be seen as making the following set asymptotically stable. 2 (ii) Path centres sub-formation: in formation, the paths of specification (i) satisfy desired relative spacial relations given by desired relative displacements q 1 , … , q n−1 ∈ ℝ 3 of the centres x i c ( ), where x i c − x i+1 c = q i , i = 1, … , n − 1. Refer to Figure 3. As  has a globally reachable node, by lemma 2 in [35], ker L = span 1. Using this, specification (ii) is equivalent to asymptotically stabilizing ⊤ . 3 (iii) Angular positions sub-formation: the agents attain, on the paths of specification (ii), desired relative on-path angular positions given by desired differences of i (relative on-path angular positions), as defined in (3). This can be alternatively given by desired differences of the heading angles i 1 . Let i − i+1 = i mod 2 , i = 1, … , n − 1, for desired angles 1 , … , n−1 , refer to Figure 3. Specification (iii) can be seen as asymptotically stabilizing It is suggested here that C x PFF can be of particular interest in applications of multi-agent systems. First, it is easy to imagine a situation where a team of autonomous robots is required to surround a compact region by following a common path, an instance of C x PFF. Applications include wild fire monitoring, search and rescue missions, and pollutants containment using flying agents. Also, in environmental monitoring in air or underwater, near surface planetary exploration, and military applications. Although a compact region can be contoured by a circular path, yet, the path class addressed in C x PFF, which is more general, could allow for better coverage. The path shape and orientation parameters are incorporated as control parameters which can be changed online. This helps to adapt the formation to changing situations, for example, spreading fires or pollutants, as well as different terrain inclinations. Also, of the control parameters are the relative path angular positions, and so uniform as well as non-uniform distributions are achievable. The ordering of the agents can be controlled as well, which can be useful if different agents perform different tasks, for example, monitoring, sampling and taking action. In addition, addressing the problem as stabilizing a goal set has an added advantage. Suppose a certain instance of C x PFF is achieved, since the design will render the goal set stable, the previous parameters can be changed on-line in small increments without the agents deviating much, during the transition, from previous states.
Second, although the case where each agent follows a different path, whether on the same or on different planes, is not familiar, at least to the author, yet one can imagine a situation where each agent is required to monitor a certain region while maintaining contact with the rest. The design presented would allow for completely arbitrary paths configurations by setting the control parameters. In this case, the on-path angular positions aspect is a part of controlling the agents relative positions. This can be useful in applications in the field of sampling and remote sensing, for example, in localization.

2.2
Approach and tools C x PFF is solved using a hierarchy of five steps of stabilizing nested sets Γ 1 ⊃ Γ 2 ⊃ ⋯ ⊃ Γ 5 , where Γ 5 equals the goal set Γ. For i ∈ {1, … , 5}, step i of the hierarchy involves designing (distributed) feedbacks to asymptotically stabilize Γ i relative to . This notion is defined as follows.
The following result, adopted from [5] gives conditions for when the previous implies that Γ is asymptotically stable. Corollary 1. For a dynamical system Σ :ẋ = f (x) with state space , a domain in ℝ n . Let l > 1 be some constant, and  ⊃ Γ 1 ⊃ ⋯ ⊃ Γ l , with Γ l non-compact, closed positively invariant sets for Σ. If for i = 1, … , l • Γ i is asymptotically stable relative to Γ i−1 , and • for some k ∈ {1, … , l }, Σ is locally uniformly bounded (LUB) 4 near k , then k , … , l are asymptotically stable for Σ. The first condition is also necessary.
Showing the second condition can be a technical challenge, refer to [3]. The following result shows that this simplifies in the case where each step of the hierarchy is solved for exponential stability. In such case, the property 'migrates' up the hierarchy, from the restriction on k−1 to the state space.
Corollary 2. For a dynamical system Σ, and Γ 1 ⊃ ⋯ ⊃ Γ l as in Corollary 1. Suppose that there exists a diffeomorphism x ↦ z such that Σ can be written aṡ In addition, Γ j , j ∈ {1, … , k} is asymptotically stable.
The proof follows from developments as in Lemma 9.4 in [36].

SOLUTION HIERARCHY
C x PFF is broken down into five hierarchic steps H1-H5.

H1-Feedback-linearizing-thrusts orientation stabilization:
Suppose (1) is fully-actuated, that is, f i 3 can take any value. In this case, the feedback linearizing thrusts with k i f ( ) > 0 and u i ( ) some C 1 design functions, render the set { ∈  :ẋ i = u i ( )b i 1 , i ∈ n} asymptotically stable, that is, the agents linear velocities are stabilized to vectors with amplitudes u i ( ) and directions b i 1 . In this step, it is required to stabilize the orientations of the rigid bodies to a direction making f i 3 = 0 for the feedback linearizing thrusts (7). This can be achieved by stabilizing the set where H2-Linear velocity stabilization: Design the feedbacks f i 1 ( ), f i 2 ( ), i ∈ n, to obtain desired linear velocities on Γ 1 , that is, stabilize the set with u i ( ) as in H1. Note that the dynamics on Γ 2 are that of steered kinematic particles.

H3-Vertical formation stabilization:
For the kinematic motion on Γ 2 , make the agents approach planes orthogonal to p 3 , with separations consistent with the desired formation. This corresponds to stabilizing the set where H4-Horizontal formation stabilization: On the planes of H3, make the agents approach the formation in specification (ii) of C x PFF. This corresponds to stabilizing the set { ∈ Γ 3 : L (3) (x c ( ) − ) = 0}. Consider the centres of rotation in (5). Define their projections orthogonal to p 3 by . Also, define the projection of orthogonal Using this, H4 requires stabilizing the set

H5-Angular positions stabilization:
While in the formation of H4, make the agents acquire the desired relative on-path angular positions, with the required speed specification. This corresponds to stabilizing This hierarchy provides extra problem features which can be of interest. As the sets expressing the different sub-problems are nested, and are to be rendered stable by control design, variations in certain formation aspects can be made while maintaining others. For example, as the design will render the agents path invariant, the formation on the path(s) can be changed without the agents leaving it. For instance, for a single path formation, the agents could be made to come close together for communication purposes, or to focus on certain regions. Also, if an agent is lost from a formation, the rest can be made to redistribute without leaving the path. Another useful property is that the design renders the path planes invariant, so, the agents can be made to follow different formations, for example, a path with different shape, individual paths on the same plane, or different Finally, the underactuation is handled at the first step of the hierarchy where on Γ 2 the system is seen as fully actuated. This implies that the specifications of the hierarchy can be changed without the need to revisit this step. Also, it is worth noting that on Γ 3 the dynamics are that of kinematic unicycles, and on account of the hierarchy available solutions for different problem for that model can be implemented, more or less directly.

CONTROL DESIGN
In this section, the solution of the previous five-step hierarchy is solved in five consecutive steps. Refer to Figure 4 for an overview.

Solution of H1
The solution of this step uses the functions u i ( ), w i 2 ( ), w i 3 ( ) to be designed in Sections 4.3, 4.4, and 4.5. Those forms will be used here directly, but on account of the hierarchy one can follow the developments in those sections first. Consider with u i ( ), w i 2 ( ), w i 3 ( ) as defined in (19), (21), (22), and  (14) and (15) render Γ 1 asymptotically stable.

Solution of H2
Consider the diffeomorphism and the candidate Lyapunov functions With (14), the time derivative of W i along the system dynamics on Γ 2 takes the formẆ i = −2k i f W i . This implies that the set Γ 2 is asymptotically stable relative to Γ 1 .

Solution of H4
The objective of this step is to design the i ( ) of (21) such that Γ 4 in (12) is asymptotically stable for (20), and i ( ) = 0 on Γ 4 . From the solution of H3, on Γ 3 . So, on Γ 4 , when i ( ) = 0, each agent follows a stationary path as required. Therefore, Γ 4 is invariant foṙ with ∈ (ℝ). In Appendix A.4, it is shown that this choice solves for the required asymptotic stability.

SOLUTION OF C x PFF
This section provides the solution of C x PFF for (1), in addition to two corollaries for fully-actuation, and degree-two underactuation with f i 1 = 0.

Degree-one underactuation
The next result follows from the design in Section 4. (1) under Assumptions 1 and 2. For i ∈ n and sufficiently small k,k,k, the feedbacks (14) where u i ( ), w i 1 ( ), w i 2 ( ) and w i 3 ( ) are given in (19), (A.6), (21) and (22), respectively, with i ( ), i ( ) andû i (b 1 ) as in (23), (25) and (28), solve C x PFF, rendering Γ 5 in (13) asymptotically stable for the closedloop system. In addition, the feedbacks render Γ 3 in (11), and Γ 4 in (12) asymptotically stable with the added properties that for initial conditions in a neighbourhood of Γ 3 the agents approach fixed planes, and for initial conditions in a neighbourhood of Γ 4 the agents converge to stationary formations (stationary path centres).

Theorem 1. Consider the n rigid bodies in
A proof sketch is provided in Appendix A.5.
Remarks: 1. The feedbacks in the theorem satisfy the information restrictions of Assumption 2. These require the knowledge of relative velocities which can be challenging to obtain in applications. On account of the hierarchy, any observer design could be incorporated, in an extra step, without modifying the original design as long as the observer errors allows for the LUB property. In this case, the actual values in the previous feedbacks are directly replaced by their estimated ones.
2. One could arrive at the previous solution by applying the backstepping algorithm in [5], developed using hierarchic set stabilization. However, using the latter directly helps to establish a relation between the technical development and solution features in regard to the formation problem.
3. As mentioned previously, the author is not aware of a work addressing this exact problem. A non-hierarchic approach would in some way involve addressing the dynamics of the variables (A.25). In this case, one would face at least two problems: how to deal with the large order of the system, and the underactuation.
4. An advantage of the hierarchic approach is that individual sub-problems could be of independent interest/use. For instance, the solution of H5 utilized that of a sub-problem in [3]. Also, some of the ideas developed in [4] has been used here as building blocks of the solution. Note that the dynamics in the solution of H3 are that of steered kinematic particles, the system of interest in [4], and that of H4 are of unicycles, the system of interest in [3]. This shows how the approach allows for previous modules to be used for solving new/more complicated problems.

Full-actuation
The next result follows directly from Theorem 1.

Corollary 3.
Under the assumptions of Theorem 1, the feedbacks (7), with constant k i f , where the functions u i ( ), w i 2 ( ) and w i 3 ( ) are given in (19), (21), and (22), respectively, with i ( ), i ( ), andû i (b 1 ) as in (23), (25) and (28), and w i 1 ( ) any C 1 uniformly bounded functions. solve C x PFF for fully actuated agents, rendering Γ ′ asymptotically stable, Remark 3. The requirement p 3 ⋅ b i 2 = 0 in Γ is not part of the problem description but rather is due to the underactuation. When the system is fully-actuated, that requirement is not necessary which renders w i 1 a degree of freedom. Due to the hierarchic nature of the approach, w i 1 can be stabilized separately to desired w i 1 ( ).

Degree-two underactuation
Suppose that for system (1), f i 1 = f i 3 = 0. In this case, C x PFF can be solved only if the desired forward speed is constant. With constant speed, specification (iii) (desired relative on-path positions) cannot be satisfied unless the path is circular. Thus, for this case C x PFF can be solved if I. the desired path is circular,ū( i ) = v (v > 0 desired speed), and ‖ ′ ( i 1 )‖ = r (r > 0 desired path radius); II. or for a convex path with u i ( ) = v, and without specification (iii) (desired on-path angular positions).
Note that in the following development u i in (21) is set equal to v in only the first equation for case I, and in the first and second for case II. In those cases, the problem can be solved using four hierarchical steps. Steps 1-3 are H1-H3, and step 4 is, for case I H4 and H5 solved simultaneously, and only H4 for case II. The design for steps 1-3 is the same as in Sections 4.1-4.3, however, showing the result for steps 2-3 is different. Consider with k i f ( ), w i 1 ( ), w i 2 ( ), w i 3 ( ) as in Theorem 1, the diffeomorphism and the functions W i in (17). The time derivative of W i along the system dynamics on Γ 1 takes the forṁ . By Barbalat's lemma [36], this implies that e i 2 (t ), e i 3 (t ) → 0. Using the facts: W i is positive definite, e i 1 = w i 3 ( )e i 2 − w i 2 ( )e i 3 , and w i 2 ( ), w i 3 ( ) are uniformly bounded, it follows that e i 1 (t ) → c i 1 where c i 1 are some constants. SinceẆ i is positive semi-definite, then for initial conditions on { ∈ Γ 1 : W i < W 0 , i ∈ n}, with W 0 sufficiently small, c i 1 < v. By this, the set Γ ′ 2 := { ∈ Γ 1 :ẋ i = v ′ i b i 1 , i ∈ n}, for some constants v ′ i > 0, is asymptotically stable relative to Γ 1 . On Γ ′ 2 , the dynamics is (18) with u i ( ) = v ′ i . Using this, the development in Section 4.2 directly applies showing that Γ ′  Remark 4. While this result is local as that of Theorem 1. Simulations have shown that the latter is more robust to large initial conditions. Still, the result might be useful in situations when applying the degree-one underactuated feedback, and the second degree of actuation is lost close to the formation.

SIMULATIONS
Consider five rigid bodies with  having a cyclic Laplacian

Formation 1
The path in Figure 5 is an ellipse with major axis of length 2, and minor axis of length 1. The path orientation is such that  = , where the major axis points in the p 1 direction.
The formation is given by (the separation between the path centres), and 1 = ⋯ = 4 = 0 (the on-path relative angular positions). Therefore, in this formation the agents should follow individual paths, neither coplanar nor concentric, belonging to horizontal planes, that is, parallel to the 1 − 1 plane, and have synchronized headings. The stabilizing feedbacks used in this simulation are those of degree-one underactuation (Theorem 1). Note that while the agents follow their respective paths, they are oriented such that their forward directions b i 1 's (marked with blue) are tangent to the paths and pointing in the direction of motion, and the unactuated directions b i 3 's (marked with red) are orthogonal to the path planes.

Formation 2
The path in Figure 6 is an ellipse similar to that of Figure 5.  The formation is such that the centres of the agents' paths form a uniform pentagon in the p 1 − p 2 plane. This is given by This gives co-planar paths on a plane that is a −30 • tilt of the horizontal around 1 . In this case too, it is required to synchronize the agents' headings, that is, 1 = ⋯ = 4 = 0. The stabilizing feedbacks used in this simulation are those of degree-one underactuation (Theorem 1). Same as in Formation 1, while the agents follow their respective paths they are oriented such that their forward directions b i 1 's (marked with blue) are tangent to the paths, pointing in the direction of motion, and the unactuated directions b i 3 's (marked with red) are orthogonal to the path planes.

Formation 3
The simulations presented in Figures 7, 8, and 9 are for a circular formation with radius r = 2, and orientation given by The agents are required to follow a common path, that is, q 1 = ⋯ = q 4 = [0 0 0] ⊤ , and the formation required is uniform, that is, 1 The simulation in Figure 7 is for degree-one underactuated feedbacks. Notice that the agents' orientations are stabilized to forward directions tangent to the path in the direction of motion, and unactuated directions orthogonal to the path plane as in the previous formations. The simulation in Figure 8 starts as in Figure 7, then after 20 s f i 1 is set to zero. The simulation in Figure 9 is for fully-actuated feedbacks. Notice that in this case the unactuated direction b i 3 (t ) does not approach the p 3 direction, but instead stabilizes at arbitrary directions. Here, w i 1 is set to zero.

PRACTICAL AGENTS
As mentioned previously in the introduction, although the agents model used here is more relevant to actual aerospace or underwater agent models, yet in order to use it in actual applications extra aspects need to be addressed such as force fields and coupled degrees of actuation. It is conjectured here that the hierarchic approach could facilitate this task. One way is as following. Instead of carrying out a new design, it is sufficient to design those agents actuator inputs to provide desired forces where d i ( ) :  → ℝ 3 is a term added to address the force fields, and angular velocities w i ( ), or torques used to achieve them. This can be done in two hierarchic steps: 1. Design actuator inputs to achieve desired actuator outputs. 2. Design actuator outputs to achieve desired forces. Refer to Figure 10. The second step involves addressing the couplings between the actuation degrees. Figure 11 shows two examples of practical agents that can be partially modelled by (1). The first is fixed wing UAV. This system usually has one force from a propeller or a different kind engine, and another force generated by the movement of the air- foil. The second system is an underactuated satellite from [29]. While satellite thrusts can be more or less than two, usually two orthogonal thrusts are needed for orbital maneuvers. One for orbital plane stabilization, and the second for stabilizing desired orbit on that plane. The work in [29] was motivated by this where only two thrusts were used to design continuous feedbacks for formation stabilization (keeping). If the satellite is already on the desired orbital plane, with the correct orientation, then only one thrust is used.
Regarding the force fields, both of the agents types above would experience gravitational force, however the way to deal with it near the earth surface, that is for a UAV, is different than that for a spacecraft. For the former, this force needs to be cancelled by the vehicle actuators. For the latter, this force cannot be cancelled, the controller should taking into account that this force shapes the dynamics of the vehicles. In [29], circular formations were addressed, and despite that it is not feasible to use the approach discussed here, that is, the two hierarchies of Figures 4 and 10, yet it was still possible to design another hierarchy to solve the problem. In addition, to gravity, UAVs would experience air resistance. To apply the previous approach, one needs to determine the function d i ( ) which is in this case the weight of the vehicle plus an air drag function. This can be challenging to specify as air drag depends on the shape and orientation of the vehicle, as well as other factors.
Regarding actuator inputs design, the first step above is usually very simple. It amounts to designing the input of an electric motor, a pneumatic valve, or a thruster … etc. to achieve a desired actuator output. The second step is the main challenge for agents like the fixed wing UAV. In that system, the second force results from the forward motion, the shape, including the angles of the control surfaces, and orientation of the vehicles. Also, the angular velocity is a function of the forward force and the angles of the control surfaces. For the satellite, the situation is simpler because usually the thrusters are fixed in certain body directions, and there are several methods for achieving desired angular velocities.

Example
Considered a quadrotor as shown in Figure 12, in earth force field, and consider the model used in [37]

FIGURE 13
Quadrotor-circular path wheref ∈ ℝ is the total thrust from the rotors, g is the gravity constant, ∈ ℝ 3 is the applied torque in body frame, parametrized by B, and J is the inertia matrix with respect to B. Suppose that a certain feedback f ( ) has been designed to achieve a certain behaviour as in the previous development. This can be used for (31) by applying the following: • Design to achieve orientation such that b 3 points in the direction of f ( ) + mg 3 . • Setf = ‖ f ( ) + mg 3 ‖. Figure 13 shows a simulation for stabilizing (31) to a circular orbit centred at [1 1 1], with its plane orthogonal to [1 0 1], and with radius value of 1. In this case, it is impossible to have, while on the path, the forward direction b 1 pointing in the direction of linear velocity, unless the path is horizontal. In the simulation, b 1 is stabilized to the path plane. Note also that in order for this design to work f ( ) + mg 3 should be well defined at least in a neighbourhood of the desired goal set.
Remark 5. Model (31) is a simplified model for a quadrotor flying vehicle. One point of simplification is the decoupled actuation degrees. Generally, the torques depend on the rotor thrusts, however, the two can be easily decoupled [38].

CONCLUSION
The paper presented a hierarchy solving a formation problem for underactuated agents in ℝ 3 , which generalizes circular formations. The problem provides certain advantages for applications of multi-agents, and the hierarchic set stabilizing approach used complements these advantages. The agents model used is close to how several practical agents can be modelled. Several questions need to be addressed to facilitate the implementation of the results in practical applications. For instance, how force fields, actuator dynamics/couplings can be handled, and how to find/extend the region of attraction of the local result. If p 1 , … , p n−1 are real, then the result follows for any k > 0. Now, suppose that some eigenvalue of L, p i , is a complex number. By the properties of L, and the Gershgorin circle theorem, p i lies in a disc, in the complex plane, centred at some r i ∈ ℝ + , with radius r i . r i is a diagonal entry of L. The roots of the term By choosing k <
These forms are equivalent to the ones in (23), (25) and (28), respectively, but are more compact. In the latter forms, the agents sensed information are expressed by variables in its body frame. Using these functions, w i 2 ( ) can be rewritten as w i provided that this angle is well defined. To this end, consider the candidate Lyapunov function Note that by the choice (15), this angle is well defined. Since b i 2 , b i 3 do not appear in i , it follows thaṫi is not a function of w i 1 . Taking the time derivative of i ,̄i along (1) giveṡ It is straight forward to show that⏞ with some constant k i w 1 > 0, the time derivative of Θ in (A.5) along (1) becomeṡ By this, and the fact that i , i are uniformly bounded, see (A.4), it follows that for initial conditions in ⊤ẋ i (t ) → 0 exponentially, and so, Γ 1 is asymptotically stable if w i 1 ( ) is well defined. By (15), this follows if u i ( )| i || i | > 0. To this end, first, k,k > 0 can be chosen small enough such that u i ( ), Let k < ≤̌( i 1 ) for all ∈ . Using this, the fact that (⋅) < 1, and the previousV i , it follows that  = { ∈  : , i ∈ n} is positively invariant. On this set, | i | > 0, and the result follows.
A.3 H3 details Consider the n × n orthonormal matrix , (A.9) and let Using (21), (22) and the fact that the dynamicṡy,̇takes the forṁ (̌( 1  1 ), … ,̌( n 1 )), and U (b 1 ) = diag(û 1 (b 1 ), … ,û n (b 1 )). Consider a change of coordinates ↦̃,̃= P −1 S (b 1 )P . By the assumption oň( i 1 ), refer to (21), this transformation is well defined for all time. Substituting this in (A.11) giveṡ Note that by (21) and (22), w i 2 ( ) and w i 3 ( ) are bounded, and hence so isṠ −1 (b 1 ). Now, consider the partition y = [ȳỹ ⊤ ] ⊤ whereȳ ∈ ℝ,ỹ ∈ ℝ n−1 . Using this in (A.12) giveṡȳ Note that since L is balanced, L ⊤ is a Laplacian of a digraph which has a globally reachable node as well. From this, and the fact that the first column of P lies in ker L, P −1 LP = diag(0, M ), for some matrix M ∈ ℝ n−1 × ℝ n−1 whose eigenvalues have strictly positive real parts, and so LPy = P[0ỹ ⊤ M ] ⊤ . Since P[ȳ 0 1×n−1 ] ⊤ ∈ ker L, Γ 3 can be rewritten as { ∈ Γ 2 : y = 0,̃= 0}, and so, Γ 3 is asymptotically stable for (18) if the origin of the (ỹ, ) sub-system of (A.13) is so. This sub-system can be viewed as a perturbed system with perturbation p( ) = for some c > 0, then, by Lemma 9.1 in [36],k could be chosen small enough such that the origin of the (ỹ, z ) sub-system of (A.13) is exponentially stable if that of the nominal system is so. The linearization of that nominal system around (ỹ,̃) = 0 takes the forṁỹ By the boundedness of S (b 1 ) andṠ −1 (b 1 ), it can be shown that the origin of this system is asymptotically stable if that is so, andk is large enough. By Lemma 1, the eigenvalues of have negative real parts except for one at the origin, and the result follows. Thus, the origin of the (ỹ,̃) sub-system of (A.13) is asymptotically (in fact exponentially) stable, for sufficiently smallk. Note that, from (A.13), and the uniform boundedness of A( ), U (b 1 ), S (b 1 ), if̃(t ) (and so (t )) converges to zero exponentially then,ȳ(t ) has a finite limit. Thus, for all initial conditions in a neighbourhood of Γ 3 , the solutions y(t ), (t ) converge to points on Γ 3 , i.e. for all such initial conditions the agents approach fixed planes.

A.4 H4 details
The following development is restricted to Γ 3 .