{'fraction_non_alphanumeric': 0.04161182001343183, 'fraction_numerical': 0.014694425789120214, 'mean_word_length': 4.1502490315439955, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 24, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multi-agent systems for quadcopters

Unmanned Aerial Vehicles (UAVs) have been increasingly used in the context of remote sensing missions such as target search and tracking, mapping, or surveillance monitoring. In the first part of our paper we consider agent dynamics, network topologies, and collective behaviors. The objective is to enable multiple UAVs to collaborate toward a common goal, as one would find in a remote sensing setting. An agreement protocol is carried out by the multi-agents using local information, and without external user input. The second part of the paper focuses on the equations of motion for a specific type of UAV, the quadcopter, and expresses them as an affine nonlinear control system. Finally, we illustrate our work with a simulation of an agreement protocol for dynamically sound quadcopters augmenting the particle graph theoretic approach with orientation and a proper dynamics for quadcopters.


1.
Introduction. Unmanned Aerial Vehicles (UAVs) present unique opportunities to investigate unsafe or otherwise unreachable locations. UAVs are already being used in numerous agricultural applications ( [20,17]), search and rescue ( [18,33]), and conservation biology ( [34]). While there are still many open issues in the use of UAVs, they also have much to offer when it comes to disaster preparation and response ( [19]).
In Spring 2018, motivated by the 2018 Kilauea volcano eruption, we conducted work in collaboration with NASA Jet Propulsion Laboratory (JPL) to develop a real-time data relay system for UAVs. A relay system augmented with an interactive visual aid for data viewing was developed and implemented to exploit the live stream, including image superposition between sequential flights to rapidly provide information about changes to the situation under study. In September 2018, we demonstrated the concept at the University of Hawai'i at Mānoa with live streaming of images from Bellows Air Force Station and video from Fissure 8 on Hawai'i Island. We then developed on-board capability for image analysis for real time situation assessment ( [12]). In parallel, we also started to study motion planning for a swarm of UAVs ( [15]), which is the problem we focus on in this paper. See [2] for a video demonstrating some of the work.
Implementing a coordinated group of UAVs, rather than a single vehicle, is naturally more efficient for many mission types. These include tracking, surveillance, and mapping ( [14]). In this regard, there are numerous reasons, driven by the variety of applications, why collective control of a group of UAVs is important. Common objectives include synchronization, formation, flocking, consensus, and rendezvous. Rendezvous is the primary focus here ( [21]). Moreover, autonomy of individual UAVs provides a flexibility and adaptability that can be beneficial in many kinds of applications. Optimization techniques, as an approach to collective control, have the general drawback of scaling poorly with the number of vehicles. We shall focus on a common alternative: the so-called graph theoretic methods. The relative simplicity of these methods makes them an attractive choice in terms of both scalability and the desire for real-time control ( [25,27]). Moreover, they apply quite naturally to situations where there is an existing network of communications between the UAVs, although the existence of such is also an additional assumption being made. The benefits include control that is decentralized in that it occurs over the network without any advanced planning, and certain objectives, such as rendezvous, can be carried out with only relative information being shared across the network.
In this paper we focus on path-planning of the most common type of UAV, quadcopters, and we focus on the dynamics of this specific type of UAV in Sections 3 and 4. Section 2, however, applies to any multi-agent UAV path configuration since we treat them as particle points. We also occasionally refer to any UAV by the common, though imprecise, more general term "drone". These small versatile vehicles have proved to be a major asset in the increasing demand for complex autonomous robotics missions, particularly those involving survey and object detection. Our end-goal is developing, for a swarm of quadcopters, the ability to autonomously make decisions based on shared information.
The outline of the paper is as follows. In Section 2 we discuss the multi-agent autonomous agreement. Agents are considered as particles and the methodology is based on graph theory. Section 3 introduces the dynamics of quadcopters and their navigation. In Section 4 we approximate our particle multi-agent trajectories by implementable ones satisfying the quadcopter's dynamic equations of motion.
It is a genuine pleasure to have our article appear in this special issue that honors Professor Bloch. The depth and breadth of his prodigious body of work inspires and enriches each of us. On those occasions when the authors feel like explorers venturing into uncharted territory, we often find a flag marking Tony's prior visitation. Indeed, Tony's contributions to the subject of geometric control of mechanical systems are keenly felt in this paper. His recent work on multi-agent formations with obstacle and inter-agent avoidance in [13] aims toward a similar goal as ours. Although the mathematical approach taken differs from the one here, some of his work in [29] does apply a graph theoretic methodology. In addition, the rigid body dynamics that we highlight in the second half of the paper is core to much of his research ( [3,5,4,30]). With sincere appreciation and gratitude, the authors offer Professor Bloch their very best wishes.

2.
Multi-Agent Autonomous Agreement. Suppose a mission with multiple agents requires that at some point the vehicles come together at a single location, for instance at the end of a mission. With batteries possibly running low, a main objective is to conserve power and rely on short-range communication. Autonomous agreement is one of the fundamental problems in multi-agent coordination, where a collection of agents are to agree on a joint state value.
2.1. Agreement Protocol. Our multi-agent systems will be viewed as networks, as described in this section. Our primary goal is to highlight the existence of an intricate relationship between the convergence properties of the so-called agreement protocol on one hand, and the structure of the underlying agent network connections on the other. Definition 2.1. We define a network of agents to be a finite collection of agents possessing some inter-connectivity via relative information exchange links. The network is modeled by a graph, denoted G(V, E), where V is the set of vertices and E is the set of edges. The agents correspond to the vertices of G(V, E) and the relative information exchange links correspond to the edges.
In the sequel we denote by V = {v 1 , v 2 , . . . , v n } the vertices corresponding to the n agents, labeled 1, 2, . . . , n. The set of edges is a subset We shall enhance this initial model to accommodate time-varying weighted networks, motivated by the desire to account for the strength of communications between agents, or possibly to take into account how the strength of communications evolves in the rendezvous problem as distances between agents change. Interestingly, suitable assumptions on the weights imply the agreement value is an invariant, and only the trajectories leading to it differ.
The agreement protocol assumes the rate of change of each agent's state is governed by a sum of relative states with respect to a subset of neighboring agents. Assume the scalar function q i (t) denotes the state of agent i for i = 1, . . . , n, and define q(t) = (q 1 (t), · · · , q n (t)) T . It follows thaṫ where N (i) is the neighborhood of a vertex v i , i.e., the collection of adjacent vertices If we adopt the notion that adjacency is symmetric, then we say the network is undirected and we can encapsulate the above system in a single matrix equation: where L(G) ∈ R n×n is called the graph Laplacian. The graph Laplacian is uniquely determined by the structure of G(V, E).
Definition 2.2. We refer to Equation (2) as the agreement dynamics, the agreement protocol, or the consensus protocol.

2.2.
Consensus Joint Value. The graph Laplacian has desirable properties when G(V, E) is a connected, undirected network. In the rest of this section we will always assume G(V, E) is a connected, undirected graph unless specified otherwise.
Proof. This is a standard result and can be found in [22,27,28].
Let Λ = diag(0, λ 2 , . . . , λ n ) and let U = [û 1û2 . . .û n ] (the hat notation indicates a unit vector) be an n × n matrix consisting of orthonormal eigenvectors of L(G) corresponding to the eigenvalue ordering above, i.e., U T U = I and L(G) U = U Λ. Since we can factor the Laplacian as L(G) = U ΛU T , the solution to Equation (2) with initial condition q(0) = q 0 is given by which can be written as: Definition 2.3. Given a set of n multi-agents, its agreement set A ⊂ R n is defined as the subspace span{ 1}. A consensus joint state value is an element of A, and the common value of its (identical) components is called its agreed state value.
The next proposition shows that when the state of a network of multi-agents satisfies the agreement protocol, the state converges to the consensus joint state value whose agreed state value is obtained by averaging the state's initial components.
Proposition 2. Let G(V, E) be a network of n agents whose state q(t) satisfies the agreement protocolq = −L(G)q. Then lim t→∞ q(t) = α 1, where α = 1 T q(0) n ∈ R is the average of the initial states. In addition, since λ 2 is L(G)'s smallest nonzero eigenvalue, its magnitude will determine the rate of convergence to α 1 ( [28,22]).
The point lim t→∞ q(t) = 1 T q 0 n 1 = α 1 ∈ R n is precisely the orthogonal projection of q 0 onto the agreement subspace and therefore minimizes the quantity q − q 0 over all possible q ∈ A, where · is the standard Euclidean distance. Furthermore, is a constant of motion since d dt for arbitrary values of t. Hence the agreed state value α = 1 T q0 n is completely determined at t = 0.
If the state of each agent is a vector quantity rather than a scalar value, the agreement protocols for the components can still be written in a single matrix equation. Indeed, assume we desire agreement on r scalar values. This means that to each agent i we associate a vector value q i (t) = (q i 1 (t), · · · , q i r (t)) and we let Q be the n × r matrix whose rows are given by the q i , and we can write the agreement dynamics as: We introduce q i (0) = (q i 1 (0), · · · , q i r (0)) the initial row vector for the state of agent i and, more importantly, q k (0) = (q 1 k (0), · · · , q n k (0)) T the initial column vector for a given state value over all agents. Since L(G) = U ΛU T , we have that the solution of (7) is given by: , which can be written as: where 1 T q k (0) n = α k ∈ R is the agreed state value for each scalar state, 1 ≤ k ≤ r. Note that here we assume the states q k , measured by the agents, are shared over the same network G(V, E). Hence the states share a common graph Laplacian, L(G). In more general situations each state variable could have a distinct network, leading to an agreement protocol of the form: 2.3. Example: Consider a collection of n agents required to meet at a single location, not given in advance, and the agents do not have access to their global positions. Rather, all they can measure is their relative distance with respect to their neighbors and they have to agree on their spatial coordinates x, y and z. In this case, r = 3 and q i (t) = (q i as an example of the agreement protocol in action. By executing the agreement protocol, the convergence to a unique rendezvous point (α x , α y , α z ) is guaranteed for any arbitrary initial conditions: (x 1 (0), y 1 (0), z 1 (0)), . . . , (x n (0), y n (0), z n (0)) (i.e., any initial spatial configuration for the n agents) so long as G(V, E) is a connected graph.

2.4.
Multi-agent Trajectories. In this section,ê j denotes the standard unit vector of all 0's with a single 1 in the jth component. We use u 1 to denote the vector 1, which is an eigenvector of L(G) with eigenvalue λ 1 = 0. With this notation, Lemma 2.4. Assume agent j is connected to every other agent, and set u 1 = 1.
Then, u 1 − nê j is an eigenvector for L(G) with corresponding eigenvalue λ n = n: Definition 2.5. We say that vertex v i is fully connected if there exists an edge between v i and every other vertex v j ∈ V, j = i. Therefore a fully connected vertex has n − 1 edges.
Recall thatû 1 = 1 √ n and λ 1 = 0. Furthermore, it is proven in [1] that the eigenvalues of L(G) are less than or equal to n and the multiplicity of the eigenvalue λ = n is equal to the number of fully connected vertices.
Theorem 2.6. Suppose vertex v j is fully connected and that q(t) satisfies the agreement protocolq = −L(G)q. Then q j (t) is a straight line.
Proof. We know that n i=1 q i = nα is a constant. We also know that: This simple ODE (Newton's Law of Cooling) implies that which is a straight line.
Note that the proof is the same regardless of whether q j is a scalar or row j of an n × r state matrix Q(t), provided it satisfies the agreement protocolQ = −L(G)Q.
In this simulation we assume the communication network between the agents is fixed and unweighted. From our algorithm, the vector of agreed state values (the agreed rendezvous position) is given by: q * = (10. 25, 9.75, 29.25).
In Figure 1 we illustrate how two different communication networks produce different trajectories, even though the agreed rendezvous position is the same for each network since it depends only on the initial spatial positions of the agents. In the first scenario (solid curve), we assume the following: agent 1 is connected to agent 2; agent 2 is connected to agents 1, 3 and 4; agent 3 is connected to agents 2, 4; agent 4 is connected to agents 2, 3. This is represented by the following graph Laplacian: In the second scenario (dashed curve), we assume: agent 1 is connected to agents 3, 4; agent 2 is connected to agent 3; agent 3 is connected to agents 1, 2; agent 4 is connected to agent 1: In the first scenario agent 2 moves along a straight line since it is connected to all other agents. In Figure 2 we compare the trajectories for agent 1 for each 2.5. Weighted Network. We can imagine a scenario where it is desirable to account for the communications strength between vehicles, or for how that strength evolves with time. Interestingly, the agreed state value is an invariant, and only the trajectories leading to it are affected [28]. Note that not all extensions are immediately straightforward, however, since the explicit solution may be difficult to write down. Not all properties can be expected to carry over; for example, the linear trajectory proof relied on cancellation due to the integer eigenvalues. We now assume that the edges are weighted, i.e., edges carry a numerical weight. If the graph G w is a weighted graph, then the off-diagonal entries of the Laplacian L(G w ) will not be limited to only 0 and −1. Instead, the −1's corresponding to edges in the graph would be replaced with −w ij = −w ji to account for the weights. The diagonal entries will once again be the absolute value of the row or column sum of the non-diagonal entries. Equations (1) and (2) therefore become: The proof that the consensus joint state value is given by α = 1q0 n 1 comes from the same logic as before since L(G w ) is a symmetric positive semidefinite matrix with eigenvalues 0 = λ 1 < λ 2 ≤ · · · ≤ λ n and L(G w ) 1 = 0. This is a consequence of the fact that weighting preserves the properties that G w (V, E) is connected and undirected. Moreover, each row and column of L(G w ) still sums to zero regardless of whether or not the edges of G w (V, E) are weighted. Hence, the values of λ 2 , . . . λ n andû 2 . . .û n , will change, but we still have λ 1 = 0 andû 1 = 1 √ n . Finally, in full generality: Theorem 2.7. Consider the dynamical system: Define the Laplacian matrix L(G w )(t) by Then the system of first order equations can be rewritten as: If the rank of L(G w )(t) is n − 1 for all t then we have: Proof. By assumption, the matrix L(G w ) is real and symmetric, therefore Hermitian. Furthermore, it is diagonally dominant with non-negative diagonal entries. Therefore it is positive semi-definite. As for the unweighted case, λ 1 = 0 with corresponding eigenvector 1. Since we assume the graph to be connected throughout the motion, the rank of L is always n − 1 and the multiplicity of the zero eigenvalue is one. We have for any time t: 0 = λ 1 < λ 2 (t) ≤ · · · ≤ λ n (t) with a corresponding set of orthonormal eigenvectorsû 1 = 1 √ n 1,û 2 (t), . . . ,û n (t).
2.5.1. Rendezvous Mission, Weighted Network. Here we imagine a scenario taking into account the strength of communications between agents. As previously stated, the agreement value is an invariant, and only the trajectories leading to it will differ compared to not taking the strength of the signal into account. We assume four agents to be initially distributed as follows: In Figure 3 we illustrate how a weighted network impacts the trajectories. For both scenarios we assume: agent 1 is connected to agent 2; agent 2 is connected to agents 1, 3 and 4; agent 3 is connected to agents 2, 4; agent 4 is connected to agents 2, 3 (as in Section 2.4.1): In the second scenario (dashed curves), we assume a weighted network with weights due to initial distance: It can be observed that for the weighted network, even though agent 2 is connected to all other agents, it does not move along a straight line anymore.
In Figure 4 we compare the trajectories for agent 1 for each of the coordinates of motion. The eigenvalues for the unweighted network are given by {0, 1, 3, 4} and the eigenvalues for the weighted network are given by {0, 13.060, 42.248, 62.558}.

2.5.2.
Rendezvous Mission, Time-varying Weighted Network. The goal of this rendezvous mission is to compare a fixed network determined by initial proximity versus a time-varying weighted network, where the weights are the inter-agent distances that vary in time, and also edges are added as the preset proximity threshold is reached.
We assume four agents to be initially distributed as follows: In Figure 5 we illustrate how a weighted network impacts the trajectories. For the fixed network (solid curves) we have: agent 1 is connected to agents 3, 4; agent 2 is connected to agent 3; agent 3 is connected to agents 1, 3; agent 4 is connected to agent 1: Initially the connectivity network is the same as for the fixed network but with weights, however as the agents get closer it changes, and all agents are connected in the end. The trajectories corresponding to the time-varying network here would be easier to approximate for a real quadcopter since they are "straighter" (the sudden jump in information leads to straighter paths toward the rendezvous, however in this case, the weights prevent the perfect cancellation that would lead to true linearity.) In Figure 6 we compare the trajectories for agent 1 for each coordinate of  3. Dynamics of quadcopters. In order to successfully implement multi-agent motion planning, we need dynamically sound equations of motion governing the movement of each drone. To this end, in this section we will derive such equations under certain realistic symmetry assumptions regarding the shape of the drone body.
Utilizing the geometric formalism of [10], we present these equations in both firstorder and second-order forms, and provide a differential geometric version in terms of the Levi-Civita connection on our configuration space equipped with a natural Riemannian metric. Finally, we give the equations explicitly in coordinates. Note that the notation in this section is independent from that in the previous section.
3.1. Setup. In this subsection we provide the basic physical notation and standard derivation of the equations of motion for a rigid body in the body-fixed frame. A rigid body can model a number of vehicles in space or in a fluid, such as satellites or underwater vehicles, which can be stabilized by spinning rotors ( [5,6]). Here, we assume our drone is nothing more than a rigid body in a fluid; later we will make certain additional assumptions explicit. We closely follow the derivations in [16,31].

The configuration space for a rigid body moving in a fluid ([3]) is
Here we equip Q with coordinates (b, R) where b ∈ R 3 represents the position of the center of mass of the body in space and R ∈ R 3×3 is a rotation matrix representing the orientation of the body aligned with the principal axes of inertia. This coordinate system corresponds to some inertial frame of reference. We regularly use the notation of the hat map, the Lie algebra isomorphism First, choose a non-inertial body-fixed frame and set the following notation: v = translational velocity in body frame Ω = angular velocity in body frame p = translational momentum in inertial frame π = angular momentum in inertial frame P = translational momentum in body frame Π = angular momentum in body frame.
We seek to derive the dynamic equations of motion for the system in the body frame. Here this will take the form of a system of ODEs in the variables (b, R, v, Ω) ∈ T Q. Our kinematic equations are:ḃ = Rv (29) We transform momenta between the inertial and body frames according to: Begin by differentiating equations (31) and (32): To rewrite these equations in a single frame, we must introduce additional quantities. In the inertial frame, let f i = external forces, i = 1, . . . , k τ i = external torques, i = 1, . . . , l, so that R t τ i = total external torque in body frame.
Now in the inertial frame, the dynamics are solely due to external forces and torques: where x i is the vector from inertial frame origin to line of action of f i . We then solve (33) This gives the evolution of the momenta in the body frame, but the equations mix momenta and velocities. In order to arrive at equations in terms of solely momenta or velocities, we must find explicit formulas relating the two. This can be acheived by the Legendre transform To make this explicit, we introduce the kinetic energies where m is the mass of the body, r C G is vector from center of gravity to body frame origin, J b is the body inertia tensor, J f is the added mass inertia tensor, M f is the added mass, C f is the added cross terms, and I 3 is the 3 × 3 identity matrix. Then by differentiating the total kinetic energy (43), we find These equations transform velocities into momenta. We now assume the body has an especially nice shape, which is realistic for most drones.
Assumption 1. Assume the body has three planes of symmetry and the principle axes of inertial coincide with the body frame axes. This implies J b , J f , and M f are all diagonal and C f = 0. When we consider gravity in Section 3.2, we also assume that this symmetry extends to the mass density of the body.
Assumption 2. Assume the center of gravity of the body coincides with the origin of the body frame. This impliesr C G = 0.
With these two assumptions, equations (44) and (45) reduce to where M = mI 3 +M f is the total mass of the system and J = J b +J f is the moment of inertia of the system. Now differentiate (46) and (47) and simplify using (39), (40), (46) and (47) to obtain the derivatives of the velocities expressed solely in terms of positions and velocities, without momenta: Finally, since for drones our fluid is air, we make the following simplifying assumption.
3.2. External forces and torques. In this section we explicitly work out the external terms E F and E T for a drone in the body frame. We have four rotors producing thrust in a quadcopter, with the ith rotor possessing angular velocity ω i producing thrust t i = K r ω 2 i , where K r is the thrust coefficient ( [26]). Here we are implicitly treating the ω i as our controls.
For the external forces, we follow [32]. The force due to drag is where v j is the jth component of the linear velocity v and C D is the vector of translational drag coefficients in the body frame. The force due to gravity is where g is the gravitational constant and the third inertial frame basis vector e I 3 represents the opposite direction of gravitational attraction. Lastly, the force produced by thrust is where the third body frame basis vector e B 3 represents the direction in which each of the motors produce thrust. Thus Ostensibly, each of these forces induces a torque, represented by the last term However, for f 1 and f 2 , no torque is induced as the line of force passes through the center of mass of the body due to our symmetry assumptions. The force f 3 , however, represents the sum of four individual forces located at each of the four rotors. Each of these forces induces a torque, and their sum gives the net torque generated by the rotors: where K d is the propeller drag coefficient and d is the distance from the drone center of mass to the rotation axis of each rotor.
The external torques are computed as in [8,32]. The torque in the body frame generated by drag is where Ω j is the jth component of the angular velocity Ω and C τ is the vector of rotational drag coefficients in the body frame. The torque due to gyroscopic effects is where J r is the moment of inertia for a rotor. Assuming the rotor has nontrivial inertia about only the vertical axis, J r has only one nonzero entry which we denote byJ r , so that J r = diag(0, 0,J r ). We then have E T = τ 1 + τ 2 .
We can now give our updated equations of motion for a quadcopter: 3.3. Coordinate expression. Here we work in coordinates η = (b 1 , b 2 , b 3 , φ, θ, ψ) on Q = SE(3) representing the position and orientation of the drone in the inertial frame. The positions b 1 , b 2 , b 3 are the standard coordinates for R 3 . The angles φ, θ, ψ are Tait-Bryan angles, known in aeronautics as roll, pitch, and yaw respectively, and sometimes referred to as Euler angles (although not the proper or classical kind). The rotation angles are given intrinsically by the sequence z − y − x or extrinsically by the sequence x − y − z. Here we take (φ, θ, ψ) ∈ (−π, π) × (−π/2, π/2) × (−π, π). These coordinates on Q induce natural coordinates for the velocities (v, where C x = cos(x), S x = sin(x) and T x = tan(x). We can write R as the composition of three pure rotations as follows. Let R i,α denote rotation by α radians about the ith principal axis of inertia. Explicitly, Then R = R 3,ψ R 2,θ R 1,φ . Now R transforms linear velocity in the body frame to linear velocity in the inertial frame, and Θ transforms the body-fixed angular velocity (Ω 1 , Ω 2 , Ω 3 ) into the Euler rate vector (φ,θ,ψ). Thus our kinematic equations (54) and (55) take the formη Substituting R into equations (56) and (57) and expanding yields the following.
Lemma 3.1. In coordinates developed above, the equations of motion for a quadcopter take the forṁ 3.4. Geometric formulation. Here we will equip Q with a natural metric allowing us to characterize the equations of motion as a second order system on Q. We then lift the dynamics to a first order system on T Q. Under our assumptions, the kinetic energy (43) can be encoded as the matrix where m and J denote the mass and inertia tensor of the drone, respectively. This is a (diagonal) positive definite symmetric matrix and induces a Riemannian metric on Q, which in turn yields the Levi-Civita connection ∇. This allows us to define the acceleration of a curve in Q as ∇γγ, where a curve γ(t) = (b(t), R(t)) has velocitẏ γ = (b(t), R(t), v(t), Ω(t)). Explicitly, we compute This connection evidently contains the fictitious forces, and is independent of choice of coordinates. The geodesics for ∇ (curves with zero acceleration) represent motions of a drone subject to no external forces. Now in general, a Riemannian metric G on a manifold Q induces the "musical isomorphisms" G : T * Q → T Q and G : T Q → T * Q. Newton's equation can be expressed geometrically as F (γ) = G ∇γγ where F is some force. In our setting, this produces the geometric version of our equations of motion: where Equation (70) represents a second order system on Q. We now express this as a first order system on T Q. To that end, we invoke the vertical lift vlft w : T q Q → T w T Q where w ∈ T q Q, and the geodesic spray S : T Q → T T Q ( [23,16]). For our system, we can explicitly compute Here the components are expressed relative to the standard left-invariant basis of vector fields Then our equations of motion can be expressed by the first order systeṁ where Γ is a curve in T Q and X is the vector field on T Q given by Next, we express this system as an affine nonlinear control system as in [3]. We can rewrite equation (71) aṡ where the first summand is the drift, and the second summand contains all the controls ω i . Explicitly, we have Note that the torque due to gyroscopic effects is linear in ω i , while the force and torque due to thrust are linear in ω 2 i . Consequently we cannot write this system as an affine control system in independent controls. As is common in the literature (e.g. [26,32]) we now simplify by ignoring the gyroscopic effects in τ 2 . This simplification does not affect the simple motion planning in Section 4.1, but the gyroscopic effects are present in the simulations of Section 4.2. Consider the controls Proposition 3. Ignoring gyroscopic effects, the quadcopter can be expressed as an affine nonlinear control system as: where the drift vector field is and the control vector fields are 4. Quadcopter multi-agent dynamics.

4.1.
Single-agent motion planning. We first explore some simple motions for a single quadcopter. We use the coordinates and notation from Section 3.
where ω 2 is a free parameter.
Proof. The desired motion keeps all variables constant except Ω 3 and ψ. This imposes ω 1 = ω 3 and ω 2 = ω 4 as well aṡ for the duration of the motion. The equations of motion in Lemma 3.1 simplify considerably; in particular, we  Then and ω 4 is a free parameter.
Proof. We desire to move only along the body x-axis with no motion in the body y or z direction, requiring the body to pitch without any yaw or roll. This imposes ω 1 = ω 3 as well asḃ 3 = φ = V 2 = Ω 3 = Ω 1 = 0 for the duration of the motion.
The equations of motion in Lemma 3.1 simplify considerably; in particular, we haveθ = Ω 2 andḃ 3 = −V 1 S θ + V 3 C θ . Using ω 1 = ω 3 and settingΩ 1 = 0 forces ω 1 = ω 3 = 1 2 (ω 2 + ω 4 ). Settingb 3 = 0 and usingθ = Ω 2 giveṡ Setting this equal to the equation of motion forv 3 , we find Denote the right side of this equation by a. Finally, using t i = K r ω 2 i , we can solve for Note that in Lemma 4.2 we could alternatively treat ω 4 as the free parameter and solve for an identical expression for ω 4 in terms of ω 2 . We have similar choices for the yaw motion in Lemma 4.1.
Moreover, a nearly identical calculation gives path planning along the second body axis. Treating ω 3 as the free parameter, we can move only in the body y direction using In order to move in another direction in the body xy-plane (with no motion in the z direction), we can simply combine the motions in In the coordinates of Section 3.3, the drone begins at (0, 0, 0, 0, 0, 0) ∈ SE(3) with motors producing the necessary amount of thrust to remain at a hovering equilibrium. Over the first 12 seconds, the drone maintains orientation while rising vertically in the body (and inertial) z-direction, ending at another hover.
Over the next four seconds, the drone maintains position while rotating π/2 radians counterclockwise about the body z-axis as described in Lemma 4.1. Finally, for the last four seconds the drone executes the motion in Lemma 4.2, flying directly along the body x-axis. Figure 7 shows the path of the drone in space as well as the orientation angles over time. Figure 8 shows the controls used to produce this motion.

Multi-agent.
Here we combine most of the work completed above in order to simulate an agreement protocol for dynamically sound quadcopters. Consider three drones on a rendezvous mission with given initial conditions. We assume each drone begins hovering with all translational and angular velocities zero. This forces each initial pitch and roll to also be zero, but we choose different initial yaws. For simplicity, we assume all three drones start in the z = 0 plane (which need not correspond to the ground) but choose different initial x, y positions. This simplifying assumption does not sacrifice much generality: if the three drones started at three different nonzero altitudes, we would simply find the average initial altitude and execute simple motions bringing each drone up or down to this height. In the coordinates (b 1 , b 2 , b 3 , φ, θ, ψ) ∈ SE(3) from Section 3.3, we choose our three drones to have initial configurations (0, 0, 0, 0, 0, 0), (0, 9, 0, 0, 0, −π/4), and (15, 9, 0, 0, 0, π/2).
For this example we will assume a complete network G = K 3 , so by Theorem 2.6 we have that all trajectories are straight lines. In particular, each drone will traverse the line segment connecting their initial position in R 3 to the rendezvous position, which is simply the component-wise average of the initial conditions: (5, 6, 0). Figure  9 shows these trajectories.
There are, however, two different methods for parametrizing these line segments and hence planning the drones' flights. First, we could use the methods from Section 2, obtaining trajectories which exponentially approach the rendezvous position. In this setting, the drones' orientations are ignored completely, as are all dynamics, giving a simple yet unrealistic path.
A main objective of this project was to investigate a second method: implement the path-planning from Section 4.1 to find controls bringing the drones along the desired line segments in a manner that respects the realistic dynamics from Section 3. This requires each drone beginning its flight by executing a pure yaw rotation to orient its body x-axis toward the rendezvous position, then fly straight in this direction. See Figure 10 to see the yaw for each drone over time; note that the total flight times differ for each drone. Also note that for each drone, the roll over time is identically zero, but the pitch is non-zero in order to produce lateral displacement (compare with Figure 7).  In order to create the desired trajectories in accordance with our equations of motion, we utilize Lemmas 4.1 and 4.2. The resulting controls -the four motor speeds for each of the three drones -are shown in Figure 11, as is the total thrust for each drone. The parameters used in this simulation were taken from [7], and appear in Table 1.
Finally, in Figure 12 we compare the paths generated by the two methods above. The solid curves represent the simple but unrealistic paths derived in Section 2. The dashed curves represent the dynamically sound paths in Figures 10 and 11. Both methods trace the same curves shown in Figure 9, but differ in their parametrizations. The top plot shows the x-coordinate over time; the solid curves for drones 1 and 2 coincide as both start with an initial x-value of 0. The bottom plot shows the y-coordinate over time; the solid curves for drones 2 and 3 coincide as both start with an initial y-value of 9. Note that the solid curves approach (but do not reach) the rendezvous position much faster than the dashed curves as they are unconstrained by orientation or realistic acceleration. To move from the simulation Figure 11. The controls used to produce the motions in Figures  9 and 10. We show the four motor speeds and total thrust for each drone as functions of time.  environment to a more realistic experiment with actual drones would require a more careful consideration of the observer design as in [9].

5.
Conclusion. The primary goals of the paper were twofold. First, we introduce the agreement protocol from a graph theoretic viewpoint and apply it to a rendezvous mission for a set of particle agents without orientation. Second, we derive the equations of motion for a quadcopter to take into account its orientation and dynamics. We then apply the quadcopter dynamics to a rendezvous mission and compare the trajectories with the ones from the network multi-agents algorithm. This is a first step towards developing more complex multi-agent systems for quadcopters. A next step will be to analyze the decoupling vector fields for the quadcopter from its affine connection control system formulation and use these trajectories to approximate geometric trajectories obtained from the multi-agents network approach. We focused here on the agreement protocol, but ongoing work shows that it can be extended to formation flying for multi-agents and more complex missions, such as search and rescue with a large team of quadcopters including leaders and followers. We are also planning to run real drone flights to demonstrate the applicability of our techniques.