On-board range-based relative localization for micro air vehicles in indoor leader–follower flight

We present a range-based solution for indoor relative localization by micro air vehicles (MAVs), achieving sufficient accuracy for leader–follower flight. Moving forward from previous work, we removed the dependency on a common heading measurement by the MAVs, making the relative localization accuracy independent of magnetometer readings. We found that this restricts the relative maneuvers that guarantee observability, and also that higher accuracy range measurements are required to rectify the missing heading information, yet both disadvantages can be tackled. Our implementation uses ultra wideband, for both range measurements between MAVs and sharing their velocities, accelerations, yaw rates, and height with each other. We showcased our implementation on a total of three Parrot Bebop 2.0 MAVs and performed leader–follower flight in a real-world indoor environment. The follower MAVs were autonomous and used only on-board sensors to track the same trajectory as the leader. They could follow the leader MAV in close proximity for the entire durations of the flights.


Introduction
Swarm robotics offer to make Micro Air Vehicle (MAV) applications more robust, flexible, and scalable (S ¸ahin, 2005;Brambilla et al., 2013).These properties pertain to a group's ability to remain operable under loss of individual members and to reconfigure for different missions.Furthermore, one can imagine that, through cooperation, a swarm of MAVs could execute tasks faster than any single MAV.The envisioned applications of such multi-agent robotic systems are plentiful.Examples of interest are: cooperative surveillance and/or mapping (Saska et al., 2016;Schwager et al., 2009a;Achtelik et al., 2012), localization of areas of sensory interest (e.g. chemical plumes) (Hayes et al., 2003;Schwager et al., 2009b), the detection of forest fires (Merino et al., 2006), or search missions in hazardous environments (Beard and McLain, 2003).In order to deploy a team of MAVs for such applications, there are certain behaviors that the MAVs should be capable of, such as collision avoidance (Coppola et al., 2016;Roelofsen et al., 2015) or leader-follower/ formation flight (Vásárhelyi et al., 2014;Cheng et al., 2014;Gu et al., 2006).These tasks are accomplished by the MAVs through knowledge of the relative location of (at least) the neighboring MAVs in the group, for which several solutions can be found in literature.
Often used are external systems that provide a global reference frame within which agents can extract their own, and the other MAVs', position.One example is Motion Capture Systems (MCSs) (Schwager et al., 2009b;Mulgaonkar et al., 2015;Kushleyev et al., 2013;Michael et al., 2010;Turpin et al., 2012;Chiew et al., 2015;Hayes and Dormiani-Tabatabaei, 2002).These systems provide highly accurate location data, but only within the limited coverage provided by the system.Alternatively, Global Navigation Satellite System (GNSS) can be used to provide similar location data (Gu et al., 2006;Saska et al., 2016;Vásárhelyi et al., 2014;Quintero et al., 2013;Hauert et al., 2011).Although GNSS is widely available, it has relatively low accuracy if com-Leader (L) C o m m u n ic a te d vL Fig. 1 Leader-follower flight with 3 Parrot Bebops, equipped with UWB modules.By estimating and communicating their relative range (R) and ego-motion (v), follower 1 ( f 1 ) and follower 2 ( f 2 ) are able to localize the leader and able to follow its trajectory with a certain time delay.
pared to MCS and therefore large inter-MAV separation is required to guarantee safe flight (Nägeli et al., 2014).Furthermore, GNSS cannot reliably be used indoors due to signal attenuation (Liu et al., 2007) and can also be subject to multi-path issues in some urban environments or forests (Nguyen et al., 2016).
To increase the versatility of the solution, MAVs should thus use on-board sensors to determine the locations of neighboring MAVs.Often, vision based methods are employed, such as: onboard camera based systems (Nägeli et al., 2014;Iyer et al., 2013;Conroy et al., 2014;Roelofsen et al., 2015), or infrared sensor systems (Kriegleder et al., 2015;Stirling et al., 2012;Roberts et al., 2012).A drawback of these systems is that they have a limited field of view.This issue can be tackled by creating constructs with an array of sensors (Roberts et al., 2012) or by actively tracking neighboring agents (Nägeli et al., 2014) to keep them in the field of view.The first solution introduces a weight penalty, while the second solution severely limits freedom of motion and scalability as a consequence of the need for active tracking of neighbors.Therefore, neither solution is ideal for MAVs; a natively omni-directional sensor would be more advantageous.One such sensor is a wireless radio transceiver.Guo et al. (2017) recently implemented an Ultra Wide-Band (UWB) radio-based system for this.Range measurements are fused with displacement information from each MAV to estimate the relative location between MAVs.However, their method suggests that each MAV must keep track of their own displacement with respect to an initial launching point.If this measurement is obtained through on-board sensors (for example, by integrating velocities) then this measurement can be subject to drift over time.Alternatively, Coppola et al. (2016) demonstrated a Bluetooth based relative localization method.Rather than using displacement information, the velocities of the MAVs, the orientation, and the height were communicated between each other, and the signal strength was used as a range measurement.
Despite the promising results of range-based solutions, a drawback of the solutions by Coppola et al. (2016) and by Guo et al. (2017) is that the MAVs need knowledge of a common frame orientation.This is established by having each MAV measure their heading with respect to North, which would be typically done with magnetometers.Magnetometers are notoriously susceptible to the local disturbances in the magnetic field.In indoor environments, disturbances upwards of 80 • can occur Afzal et al. (2010).The difficulty of establishing a reliable direction towards North in an indoor environment is a well known problem.Solutions are found in the form of complementary filters (Roetenberg et al., 2005(Roetenberg et al., , 2007;;Afzal et al., 2011;Yuan et al., 2015), or the use of redundant magnetic sensors to compensate the local disturbances (Afzal et al., 2010;Li et al., 2006).These solutions, however, may be unnecessary for the purpose of relative localization, since a shared reference frame is not theoretically necessary when performing range based relative localization (Zhou and Roumeliotis, 2008;Martinelli and Siegwart, 2005).
The main contribution of this paper is an analysis of the consequences of removing the heading dependency in range based relative localization, leading to the development and implementation of a heading-independent relative localization and tracking method that is accurate enough for full onboard indoor leader-follower flight, as shown in Fig. 1.The analysis is provided by a formal observability analysis and by performing limit-case simulations.Differently from the work of Zhou and Roumeliotis (2008) and Martinelli and Siegwart (2005), the analysis also considers the inclusion of acceleration information, since this is commonly known by MAVs from their Inertial Measurement Unit (IMU).Furthermore, our analysis specifically focuses on the implications of removing a heading dependency on the performance of the relative localization filters and on the relative maneuvers that the agents can perform in order to guarantee that the filter remains observable.The observability analysis will show that the task of leader-follower flight is especially difficult with range-based relative localization methods, because it does not allow for the MAVs to fly parallel trajectories.We then use the insights gathered for the development and implementation of a heading-independent leader-follower system that we are able to use on-board of autonomous MAVs operating indoors.The MAVs rely only on on-board sensors, using UWB for both communication and relative ranging.
The structure of the paper is as follows.First, in Sect.2, we compare the theoretical observability of range based relative localization systems both with and without a reliance on a common heading.The findings from Sect. 2 are verified through simulation in Sect.3, where we also evaluate the difference in performance that can be expected.We carry this information forward in Sect.4, where a headingindependent system is implemented on real MAVs, and where we show the results of our leader-follower experiments.The results are further discussed in Sect. 5. Finally, the overall conclusions are drawn in Sect.6. Future work is discussed in Sect.7.

Observability of the Relative Localization Filter
In this section, an observability analysis is performed that specifically focuses on the practical implications of performing range based relative localization both with and without reliance on a common heading reference.The purpose of the eventual relative localization filter is for an MAV (say MAV 1) to be able to track the position of another MAV (say MAV 2).Despite our focus on MAVs in particular, the conclusions that follow hold for any general system that can provide the same sensory information.Furthermore, the results can be extrapolated to more than two MAVs, as will be demonstrated in Sect. 4. For clarity, only MAVs 1 and 2 are considered in the coming analysis.

Preliminaries
We will conduct the analysis by studying the local weak observability of the systems (Hermann and Krener, 1977).With an analytical test, briefly introduced in the following, local weak observability can be used to extract whether a specific state can be distinguished from other states in its neighborhood.
Consider a generic non-linear state-space system ∑: The system ∑ has state vector x = [x 1 , x 2 , . . .x n ] ∈ R n , an input vector u ∈ R l , and an output vector y ∈ R m .The vector function f(x, u) contains the definitions for the time derivatives of all the states in x and the vector function h(x) contains the observation equations for the system.The Lie derivatives of this system are: Where ⊗ is the Kronecker product and ∇ is the differential operator, defined as Fig. 2 Reference frames used in this paper.Frame I in purple is the earth-fixed North East Down frame (assumed to be inertial).Frames H 1 (blue) and H 2 (red) are body fixed reference frames for MAVs 1 and 2, respectively.
h.Using these definitions, an observability matrix O can be constructed, as in Eq. 6.
A system is locally weakly observable if the observability matrix is full rank (Hermann and Krener, 1977).

Reference Frames
For the analyses that follow, consider the reference frames schematically depicted in Fig. 2. Denoted by I is the Earthfixed North-East-Down (NED) reference frame, which is assumed to be an inertial frame of reference.Denoted by H i (i = 1, 2) is a body-fixed reference frame belonging to MAV i.Its origin is coincident with MAV i's centre of gravity, and its location with respect to the I frame is represented by the vector p i .H i is a horizontal frame of reference, such that the z-axis of the H i frame remains parallel to that of the I frame.The H i frame, however, is rotated with respect to the I frame about the positive z-axis by an angle ψ i , where ψ i is the heading that MAV i has with respect to North, also referred to as its yaw angle.The rate of change of ψ i is represented by r i .Note that the H i frame is different from a typical bodyfixed frame B i , which uses the three Euler angles for roll, pitch, and yaw to represent its orientation with respect to the I frame.The reason for using H i rather than B i is that it simplifies the kinematic relations without having to impose additional assumptions, such as the roll and pitch angle of the MAV being small.

Nonlinear System Description
We shall study the case where MAV 1 attempts to estimate the relative position of MAV 2. We use p to denote this relative position, such that p = p 2 − p 1 (see Fig. 2).Furthermore, let v i and a i be the linear velocities and accelerations of frame H i with respect to frame I expressed in frame H i , respectively.Finally, let ∆ ψ represent the difference in heading between MAVs 1 and 2, such that ∆ ψ = ψ 2 − ψ 1 .
From this point on, we shall assume that the MAVs are capable of measuring their own height.Since the horizontal plane of H i matches the horizontal plane of I, the height can be treated as a decoupled dimension that does not influence the observability, provided that it is measured.Therefore, for the sake of brevity, the height is not included in the system description.The vectors for the relative position p, the velocity v i , and the acceleration a i can thus be expanded as The rate of change of ∆ ψ is ∆ ψ = r 2 − r 1 .Note that the value for r i is not equal to the yaw rate as would commonly be measured by an on-board rate gyroscope in the body frame B i .Instead, r i is expressed as: where qi and ri represent the true pitch and yaw rate as would be measured by a rate gyroscope, and φ i and θ i are the roll and pitch angles of the MAV.However, for the sake of simplicity, r i will be referred to as the MAV's yaw rate.Similarly, a i , which is the value for the linear acceleration of the H i frame expressed in coordinates of the H i frame, is not equal to what is measured by the on-board accelerometer.Instead, it is equal to: where s i is the specific force measured in the body frame B i by the accelerometer of MAV i.Furthermore, c(α) and s(α) represent short hand notation for cos(α) and sin(α), respectively.The matrix in this equation consists of the first two rows of the rotation matrix from the B i frame to the H i frame.
Following the above, the complete state vector of the system is given by x = [p , ∆ ψ, v 1 , v 2 ] , and the input vector is u = [a 1 , a 2 , r 1 , r 2 ] .The continuous time state differential equations can be written as: where R is the 2D rotation matrix from frame H 2 to H 1 : The matrices S 1 and S 2 are the skew-symmetric matrix equivalent of the cross product, adapted to the 2D case.The matrix S i is equal to: The variables a i and r i are inputs into the system and MAV 1 must thus have knowledge of these values.However, these are typically available from accelerometer and gyroscope data in combination with the appropriate relations given in Eq. 7 and Eq. 8.
Finally, Eq. 9 needs to be complemented with an observation model.Apart from the height, which must be measured but is not included in this analysis, the MAVs should be able to measure the relative range between each other, along with their own and the other's velocities.Then, the analysis that follows aims to study the difference between the following two scenarios: a scenario where the above measurements are the only measurements and a scenario where the MAVs are additionally capable of observing each other's headings.The situation where the MAVs can observe a heading is referred to as ∑ A and the situation where a heading is not observed is referred to as ∑ B .
∑ A : The scenario where ψ 1 and ψ 2 are observed is equivalent to ∆ ψ (the difference in headings) being observed.Therefore, for ∑ A , the observation model is: Note that the observation equation h A1 (x) is slightly modified with regards to the previously mentioned measurements.Rather than observing the range between the two MAVs (i.e.||p|| 2 ), half the squared range is observed (i.e. 1 2 p p).This change makes the observability analysis more convenient without affecting its result.Both ||p|| 2 and 1 2 p p contain the same information as far as observability of the system is concerned (Zhou and Roumeliotis, 2008).∑ B : In this case, the headings of the MAVs are not measured, and it is thus not not possible to observe the difference in heading ∆ ψ directly.For ∑ B , the observation model is: The effect of the difference in the observation equations is studied in the following sections.

Observability Analysis with a Common Heading Reference
For system ∑ A , which uses the observation model from Eq. 12, the first entry in the observability matrix is equal to: where I nxn represents an identity matrix of size nxn and 0 mxn represents a null matrix of size mxn.We can already deduce simplifying information from Eq. 14 that will aid the subsequent analysis.First, note that, for the higher order terms in the observability matrix, the last 5 columns do not contribute to increasing its rank, because these columns are populated with an identity matrix.Furthermore, these higher order terms in the observation matrix (corresponding to the observations of ∆ ψ, v 1 , and v 2 ) only have terms in those last 5 columns because none of the higher order Lie derivatives corresponding to those observations depend on the state p.For this reason, these need not be computed and we can thus omit them for brevity.The remainder of this analysis considers only the terms corresponding to observation h A1 (x) = 1 2 p p. The first order Lie derivative corresponding to the observation h A1 (x) = 1 2 p p is equal to: Next, remembering that S 1 is a skew symmetric matrix, such that S 1 + S 1 = 0 2x2 , the following identity is obtained: Using this identity, it is can be verified that the second term in the observation matrix corresponding to h A1 (x) is: At this point, it would be possible to continue calculating higher order terms for the observability matrix, but in practice this is not necessary.The first term of the observability matrix as shown in Eq. 14 already presents a matrix of rank 6.Since the state is of size 7, this means that only 1 more linearly independent row needs to be added to the observability matrix to provide local weak observability of the system.Furthermore, it is of practical interest to study the scenarios in which the system is locally weakly observable with a minimum amount of Lie derivatives involved in the analysis.This is due to the fact that in practice all signals are noisy, and differentiation of a noisy signal inevitably leads to increasingly noisy signals.It will be demonstrated that the terms presented in Eq. 17 are sufficient, under certain conditions, to make the observability matrix full rank.
As mentioned, Eq. 14 already shows that the last five columns of the observability matrix are no longer of interest to increase its rank.Furthermore, only the observation of h A1 (x) = 1 2 p p provides non-zero terms in the first two columns of the observability matrix.Therefore, the following matrix can be constructed by collecting the terms of the first two columns in the observation matrix belonging to observation h A1 (x): where the first term is from the zeroth order Lie derivative (see Eq. 14) and the second term from the first order Lie derivative (see Eq. 17).The system is thus observable with a minimum amount of Lie derivatives if the matrix given by Eq. 18 has two linearly independent rows.By the definition of linear independence, this means that the following condition must hold to guarantee local weak observability of the system: where c is an arbitrary constant.
The condition in Eq. 19 essentially tells us that the relative velocity of the two MAVs should not be a multiple of the relative position vector between the two.For more practical insight, we can extract more intuitive conditions that must also be met for Eq.19 to hold.These conditions are: The first condition tells us that the x and y coordinates of the relative position of MAV 2 with respect to MAV 1 should not be equal to 0. In practice, this would only be possible if the MAVs were separated by height, for otherwise their physical dimension would prevent this condition from occurring.The second condition tells us that one of the two MAVs needs to be moving to render the filter observable, and that the observability is indifferent to which of the MAVs is moving (hence the or operator).The third condition tells us that the MAVs should not be moving in parallel at the same speed (note the rotation matrix R that transforms v 2 to the H 1 frame).Whilst these three conditions are easier to consider, it should be noted that they form only a subset of the conditions imposed by Eq. 19.For example, the scenario where MAV 2 is stationary, and MAV 1 flies straight towards MAV 2, does not violate any of these three conditions.It does, however, violate Eq. 19.Therefore, the observability of a state and input combination should be checked against the full condition in Eq. 19.

Observability Analysis Without a Common Heading Reference
After determining the conditions under which system ∑ A is locally weakly observable, we compare it to the system where the heading measurements are no longer present.We now consider system ∑ B , whose observation equation (Eq.13) does not include the state ∆ ψ.For this system, the first term in the observability matrix is: Eq. 23 is very similar to Eq. 14, but with the important difference that the row corresponding to the observation of ∆ ψ is null.Consequently, the matrix is only of rank 5, rather than rank 6.Since the state size is still 7, a minimum of two more independent rows must be added to the observability matrix to make the system locally weakly observable.Once again only the terms corresponding to the observation h B1 (x) = 1 2 p p have terms that could increase the rank of the observability matrix.This means that this time a minimum of two more Lie derivatives must be calculated.
It can be verified that the first derivative L 1 f h B1 , and thus its state-derivative ∇L 1 f h B1 , are exactly the same as for ∑ A .Therefore, these need not be calculated anymore and are given by Eq. 15 and Eq.17, respectively.The second order Lie derivative is equal to: Some terms in Eq. 24 can be seen to drop out when the equation is expanded.For example, the yaw rate of MAV 1 (r 1 ) cancels out completely.Therefore, Eq. 24 reduces to: The state derivative of L 2 f h B1 can then be shown to be equal to Eq. 26.Once again, note that some terms drop out (this step has been omitted for brevity).
Just as for ∑ A , a part of the observation matrix can be extracted for analysis.This time, the first three columns in the observation matrix (as opposed to two) are collected for the observation h B1 (x) = 1 2 p p. Also, this time the terms up to and including the second order Lie derivative are minimally needed to obtain a full rank observability matrix.The following matrix is obtained: In this case, obtaining the conditions for which this is a full rank matrix is less obvious due to the plethora of terms.Rather than directly demonstrating linear independence of the three rows in Eq. 27, the determinant |M B | may be computed and demonstrated to be non-zero.This is done as follows.Recall that p = [p x , p y ].Furthermore, suppose Then, matrix M B can be written as: The determinant of M B can be computed using a cofactor expansion along the last column of M B .This results in: Now, the following identity can be used: where and [p x , p y ], the determinant of M B becomes: This can be simplified and written as: This system is thus locally weakly observable with a minimum amount of Lie derivatives if |M B | is non-zero.Due to the specific properties of the A matrix in this determinant (see Eq. 30), the following equation must hold to render the determinant |M B | non-zero: where k is an arbitrary constant.
It is difficult to find an intuitive interpretation for Eq.33.Just as for Eq. 19, we can extract a more intuitive subset of conditions that also definitely must be met for the system to be observable.These conditions are: where s an arbitrary constant.
The first condition tells us that the determinant |M B | is zero if the x and y coordinates of the origins of frames H 1 and H 2 coincide.This is the same as for ∑ A .The second condition tells us that both MAVs need to be moving.This movement may be either through having a non-zero velocity, or through having a non-zero acceleration (the violation of which is shown in Fig. 3a).The third condition tells us that the MAVs may not move in parallel, as in Fig. 3b, unless at least one of the MAVs is also accelerating at the same time.Note that this time the MAVs are not allowed to move in parallel regardless of whether they are moving at the same speed or not (notice the scalar multiple s).By comparison, the equivalent condition for ∑ A only specified that the MAVs may not move in parallel at the same speed.
In order to study these intuitive conditions in further detail, we evaluated how the observability of the system is affected once the relative position p between the MAVs changes.By varying the p x and p y values of the vector p around the originally set values for p (as in Fig. 3), we analyzed the observability of the system for different relative positions, while keeping the velocities and accelerations constant.The measure for observability was obtained by interpreting the meaning of Eq. 33.It essentially tells that the left hand side of the equation should not be parallel to the relative position vector p.Therefore, a practical measure of observability is how far away the left hand side of equation Eq. 33 is from being parallel to p, which can be tested with the cross product.The absolute value of the cross product is then used as a measure of the observability of the system.This paper con- Fig. 3 Representations of four unobservable state and input combinations.The relative position p, the velocities v i , and the accelerations a i of MAVs 1 and 2 are depicted.siders a cross product less than a value of 1 to be unobservable. 1or the case of the second (Eq.( 35), Fig. 4a) and the third intuitive condition (Eq.( 36), Fig. 4c) it can be seen that a varying p does not affect the unobservability in the color map.Once an acceleration vector is added to the state of MAV 1 in both cases, specifically a 1 = [0.30.3] , the color plots show that for a set of relative positions, the system does become observable again.However, the chances of the MAVs ending up in an unobservable state are still significant within an operating area of 100 m 2 .
The three intuitive conditions we extracted are only a subset of all conditions imposed by Eq. 33.This means that there exist state and input combinations that satisfy the three intuitive conditions, but that do not satisfy Eq. 33.In order to study what the implications of the full unobservability condition in Eq. 33 are, we used the Nelder-Mead simplex method to find other points in the state and input space that violate the full observability condition.Two examples are shown in Fig. 3c and Fig. 3d.These scenarios do not violate any of the intuitive conditions given by Eq. 34-36.The relative position is non-zero, both MAVs have non-zero ve- locities and accelerations, and the velocity vectors are not parallel.Nevertheless, they violate Eq. 33.Based on this, color maps for the unobservable conditions in Fig. 3c and Fig. 3d are given in Fig. 4e and Fig. 4f, respectively.
Both color maps of Fig. 4e and Fig. 4f clearly show a non-linear relationship between the relative position vector p and the observability of the system.Moreover, both maps show a different non-linear relationship.Fig. 4e shows more of a hyperbolic relationship, whereas the unobservable region in Fig. 4f looks more elliptical.It can be shown that different conditions show yet other relationships between the observability of the system for different relative positions p.Moreover, these relationships only show what happens in two dimensions (for the two entries in the vector p).In reality, the observability condition in Eq. 33 presents an 11 dimensional problem.It is therefore still difficult to deduce general rules from these results.What the latter two color maps do have in common is that the unobservable relative positions are in all cases vastly outnumbered by the observable relative positions.This is different than what was observed for situations that would violate any of the more intuitive conditions in Eq. 35 and Eq.36.

Comparison of the Two Systems
Finally, the results from the observability analysis of both systems will be compared.This will answer the question of what practical implications there are when moving from a system that relies on a common heading reference to a system that does not.
A primary result of the analysis is that removing the relative heading measurement results in a system that requires at least one extra Lie derivative in the range observation to make the system locally weakly observable.This is an important result, because it tells us that the headingindependent system ∑ B relies more heavily on the range equation than ∑ A .Without a heading observation, the range measurement serves to estimate a total of three states, as opposed to two in ∑ A .Some of this information is contained in the second derivative of the range observation, and it is a well known fact that derivatives of a noisy signal become increasingly noisy.In practice, this means that any system that wishes to perform range-based relative localization without a heading dependency needs an accurate and low-noise range measurement.
Another important result is that the criteria posed for ∑ B specify that both MAVs must be moving.Contrarily, the criteria for ∑ A specify that only one of the MAVs must be moving.Whilst this result might not be as relevant for MAV teams (as the MAVs will typically be moving anyway), this result can be important for other applications of range-based relative localization.Think, for example, of the case where a single static beacon is used to estimate the position of a flying MAV using only range sensing and communication.The results of our analysis show that ∑ B is not observable in this case, and thus a common heading reference must be known for such a system to work (or, alternatively, the MAV must track the beacon and then communicate its estimate back to the beacon).Note that, in the case where one of the participants is not moving, it can be shown that even the higher order Lie derivatives in ∑ B will not succeed in making the observability matrix full rank, so that this statement generally holds.
A third difference is found in the condition for parallel movement of the two MAVs.∑ A requires that the MAVs should not move in parallel at the same speed (which can be translated to mean that there should be a non-zero relative velocity between the two MAVs).Instead, ∑ B requires that the MAVs should not be moving in parallel regardless of speed.Therefore, even if the second MAV were to be moving twice as fast as the first, the filter would not be observable as long as the direction of movement is the same.
However, ∑ B , can bypass this condition in some cases if either of the MAVs is also simultaneously accelerating.Similarly, it can be shown that ∑ A is able to bypass the parallel motion condition with acceleration, although a second order Lie derivative would be necessary in that case.

Verification through Simulations
In this section, we further investigate the conclusions drawn from the analytical observability analysis.At first, a kinematic, noise-free study is performed to verify and confirm the differences in the observability conditions for ∑ A and ∑ B .Afterwards, the influence of noise and disturbances on the filter are studied.

Filter Design
The filter of choice, used throughout the rest of this paper, is an Extended Kalman Filter (EKF).This choice was made because this type of filter fits intuitively with how the statespace system was described in Sect. 2. The EKF also uses a state differential model and an observation model.The state differential model can thus be kept exactly as the one given earlier in Eq. 9.The observation models for ∑ A and ∑ B are also kept almost the same as given in Eq. 12 and Eq. 13, with the only adjustment that mow the full range ||p|| 2 is observed, rather than half the squared range 1 2 p p. Furthermore, using the EKF is in line with earlier research on rangebased relative localization (Coppola et al., 2016).
An EKF has parameters that need to be tuned, namely: the initial state, the system and measurement noise matrices, and the initial state covariance matrix.The initial state is an important setting that will be described where appropriate in the next sections.The matrices are always tuned to correspond to the actual expected values.The measurement noise matrix is tuned based on the expected noises on the measurements, and similarly for the system noise matrix.However, since some of the simulations also make use of perfect measurements (with zero noise) and since a noiseless entry in the measurement noise matrix is not possible, the corresponding entries are then given a small value of 0.1.

Kinematic, noise-free study of unobservable situations
In the first simulated study, the two MAVs that are studied have kinematic trajectories that can be described analytically.The MAVs also have perfect noise-free knowledge of the inputs and measurements.The kinematic and noisefree situation is used to confirm conclusions drawn in the observability analysis performed in Sect. 2.
The two MAVs involved in the EKF are designated MAV 1 and MAV 2. MAV 1 shall be the host of the EKF and shall attempt to track the relative position of MAV 2. For clarity, this MAV is denoted as the host of the filter.MAV 2 is the one whose position is tracked by MAV 1.It does not run an EKF.For clarity, this MAV is denoted as the tracked MAV.The following three scenarios are studied: 1. MAV 1 (host) is moving and MAV 2 (tracked) is stationary.2. MAV 1 (host) is stationary and MAV 2 (tracked) is moving.3. MAV 1 (host) and MAV 2 (tracked) are both moving in parallel to each other at different speeds.
These scenarios have been chosen because they match the intuitive conditions where ∑ A is observable, but ∑ B is not.These are limit cases and therefore provide valuable verification of the analytically found differences between the two systems.
The simulations will show whether these different scenarios have convergent EKFs or not.The focus of this analysis is on the estimation of the relative position p and the relative heading ∆ ψ.Since the velocities are observed directly, these are observable regardless of the situation, and are thus not shown.
The initial velocities of MAVs 1 and 2 are initialized to their true value, since these are not the variables of interest in this analysis.The initial position and relative heading are initialized with an error, the specifics of which will be given in the respective scenarios.The yaw rates and headings of both MAVs are kept at 0 rad/s and 0 rad, respectively.The EKF runs at a frequency of 50 Hz.
The error measure throughout this paper is the Mean Absolute Error (MAE).The separate x and y errors in the relative location estimate p are combined according to the norm ||p|| 2 .This choice was made because the separate errors in x and y directions offer little additional insight and are mostly very similar.

MAV 1 (host) moving, MAV 2 (tracked) stationary
Previous analytical analysis has shown that ∑ A is locally weakly observable, while ∑ B is not observable.This result is therefore expected to be reflected in the simulation as well.
As can be seen in Fig. 5, both the relative position p error and the relative heading ∆ ψ error quickly converge Fig. 6 ∑ B EKF convergence for case 1: MAV 1 (host) moving, MAV 2 (tracked) stationary to 0. Contrarily, the observability analysis of ∑ B has shown that this scenario is not locally weakly observable, because the second condition is violated, i.e., one of the MAVs is not moving.However, Fig. 6 shows that the ||p|| 2 error converges to 0 just as rapidly as for ∑ A .A more thorough inspection shows that the unobservable state of the system is in fact ∆ ψ, which is the one that does not converge.This is a favorable result, since the relative position is typically the variable of interest, rather than the difference in heading.
The reason that this occurs lies in the information provided by the first state differential equation.This equation tells us that ṗ = −v 1 + Rv 2 − S 1 p.The only dependency that this equation has on the relative heading ∆ ψ is in the rotation matrix R. Therefore, as long as v 2 is equal to 0, the differential equation for ṗ has no dependency on the relative heading between the two MAVs.The convergence of p therefore remains unaffected.The situation changes when it is v 2 that is non-zero and v 1 that is zero.This case is studied next.

MAV 1 (host) stationary, MAV 2 (tracked) moving
For this case, all of the parameters are the same as for case 1, with the only difference being that now v 1 = 0 and v 2 = [1, 0] .The analytical observability analysis has shown that this scenario is locally weakly observable for ∑ A .As expected, it can be seen in Fig. 7 that both the errors for p and ∆ ψ converge rapidly to 0. The observability analysis has then shown that ∑ B is not locally weakly observable in this scenario.Indeed, Fig. 8 shows that both ||p|| 2 and ∆ ψ do not converge and that ||p|| 2 even diverges.This time, because v 2 is not equal to 0, the state differential equation for the relative position of MAV 2 has a dependency on the relative heading state ∆ ψ.Because ∆ ψ does not converge to its true value, and eventually settles at an error of approximately 1.5 rad, there is a large inaccuracy in the state differential equation for ṗ.This consequently results in an ever increasing error in p, since MAV 1 essentially 'thinks' that MAV 2 is flying in a different direction than it really is.
This shows the reason as to why it is generally not possible for a stationary vehicle (or beacon) to be tracking a moving vehicle using range-only measurements and velocity information without a common heading reference.Contrarily, it is possible for a moving vehicle to be tracking a stationary vehicle or beacon's position.This is entirely caused by the fact that a vehicle will always be 'aware', in its own body frame, of the direction it is moving in and hence does not need a convergent estimate of the relative heading with respect to the vehicle it is tracking.However, when the vehicle it is tracking does move, it needs this convergent estimate of the relative heading to know which direction the other is moving in.

MAV 1 (host) and MAV 2 (tracked) moving in parallel at different speeds
Finally, the case where both MAVs are moving in parallel, but at different speeds, is studied.Once more, most of the parameters are kept the same as those presented under case 1.This time, the velocity of MAV 2 is set to v 2 = [1, 0] According to the observability analysis, this is one of the limit cases where ∑ A is still just observable, but ∑ B is not.Indeed, Fig. 9 shows convergent behavior for ∑ A , whereas Fig. 10 shows divergence for ∑ B .Note that the filter for ∑ B has a decreasing error in ∆ ψ.However, the convergence for ∆ ψ is very slow (notice how this situation has been simulated for a much longer time than the previous cases).Furthermore, the error for p continues to rise indefinitely.
This result concludes the noise-free simulations that compare the performance of the filters for ∑ A and ∑ B .These simulations verify that the conclusions regarding the differences between the two filters in Sect. 2 also hold true when translated to a simulation environment.

Kinematic, noisy range measurements study of observable situation
Whilst a noise-free study demonstrates the feasibility of the proposed filter and can verify the differences between ∑ A and ∑ B , it is also important to study the filter's performance when presented with noisy data.Not only is this more representative of the filter's performance in practice, but it also can be used to verify one of the main conclusions that were drawn in the observability study, namely that ∑ B needs information present in the second derivative of the range data to be observable, compared to only a first derivative for ∑ A .It is consequently expected that, with all other parameters fixed, ∑ B will perform increasingly worse as the range data becomes more noisy.In this study, we steer away from unobservable scenarios.The intent now is to study both filter's performances for the case where the filters are known to be observable, in order to compare their performance.For this reason, the trajectories of MAV 1 (host) and MAV 2 (tracked) are designed so as to stay clear of the unobservable situations and to excite the filter properly through relative motion.The trajectories that we devised for this study are perfectly circular, and we assume that the MAVs fly at the same height.
The trajectories, depicted in Fig. 11, can be described in polar coordinates [ρ, θ ].MAV 1 flies a circular motion at an angular velocity θ1 = ω 1 with radius ρ 1 , and MAV 2 flies at angular velocity θ2 = ω 2 with radius ρ 2 .To ensure that both MAVs have sufficient relative motion, one MAV flies clockwise and the other counter clockwise, such that ω 1 = −ω 2 .Moreover, the radius of MAV 2's trajectory is 1 meter larger than MAV 1's trajectory, and is offset by 90 • in angle, such that ρ 1 = ρ 2 − 1 and θ 1 = θ 2 + π 2 .The radius difference in the trajectories ensures that the situation p = 0 is avoided, and the angle offset ensures that the relative velocities are distributed more or less equally in x and y directions.Note that, for simplicity, both MAVs keep a steady heading such that ψ 1 = ψ 2 and r 1 = r 2 = 0. Switching back to Cartesian coordinates, the trajectories can thus be analytically described as follows.MAV 2's position vector in time is given by: MAV 1's position vector in time can be described by: The equations for v i (t) and a i (t) can be obtained by taking the time derivatives with respect to p i (t), i = 1, 2. Note that this is not true for the general case, since H i is a rotating frame of reference, but in this case it is possible because the MAVs keep a constant heading equal to 0 rad.By setting ρ 2 = 4 m and ω 2 = 2π 20 rad, the trajectory of MAV 2 becomes a circle with a radius of 4 m that is traversed in 20 s.To comply with the previously defined constraints, ρ 1 and ω 1 are 3 m and − 2π 20 rad/s, respectively.These values are representative of what a real MAV should easily be capable of and result in relative velocities of about 1 m/s in x and y directions between the two MAVs.
The study will test the performance of the relative localization filter as seen from the perspective of MAV 1, who is thus tracking MAV 2. The filter is fed perfect information on all state and input values, except for the measurement of the range ||p|| 2 between the two MAVs.The range measurement are artificially distorted with increasingly heavy Gaussian white noise.The measured range fed to the filter is thus ||p|| 2,m = ||p|| 2 + n(σ R ), where n(σ R ) is a Gaussian white noise signal with zero mean and standard deviation σ R .The standard deviations that are tested are 0 (noise free), 0.1, 0.25, 0.5, 1, 2, 4, and 8 m.In practice, a standard deviation of 8 m could be consider quite high, but this is intentionally chosen with the intent to observe a significant difference in the error.Since this study keeps all the other measurements and inputs noise free, the noise on the range measurement needs to be higher to get a significant increase in the localization error.
This time the EKF runs at 20 Hz, which is more representative of our real-world set-up, discussed later in Sect. 4. The described flight trajectory is simulated for 20 seconds each run (which is thus one complete revolution of the circular trajectory).The EKF is initialized to the true state to exclude the effects of initialization.
For each particular noise standard deviation, both the filter for ∑ A and for ∑ B are simulated with 1000 different noise realizations.For each realization the MAE of the estimated p with respect to its true value is computed, again by considering the combined error in the estimate of ||p|| 2 .After 1000 realizations, the Average MAE (AMAE) is computed to extract the average performance for all noise realizations.
The resulting AMAE values for systems ∑ A and ∑ B are given in Tab. 1 and are plotted in Fig. 12.As expected, at very low noise values on the range measurement, both the filters for ∑ A and ∑ B have very similar error performance.With no noise on the range measurements, the difference between the two filters is only 4 mm.However, since the filter for ∑ B is more sensitive to noise on the range measurements, it quickly starts to perform worse than ∑ A as the noise on the range measurement is increased.
This result is in line with the analytical results presented in Sect. 2. However, it also raises the question of whether removing the dependency on a common heading reference poses any advantage, since ∑ A performs consistently better than ∑ B .The reason for this result lies in the fact that the studied scenario uses perfect measurements for all the sensors except for the measured range.As mentioned in the introduction, the heading observation is notoriously troublesome and unreliable, especially in an indoor environment (Afzal et al., 2010).Therefore, it would be valuable to study what would happen to this analysis in the case where the heading estimate is not perfect.This is presented next.

Kinematic, noisy range measurements, and heading disturbance study for observable situation
In order to compare the results obtained with an imperfect heading measurement to those obtained in the previous section, the same trajectories are simulated (as in Eq. 38 and Eq.37 for MAVs 1 and 2, respectively).All the other simulation parameters are also kept the same, with one exception.This time, a disturbance is introduced on the heading measurement.The simulated disturbance is modeled to look similar to how a real local perturbation in the magnetic field would perturb a heading estimate.The actual magnetic perturbation and the corresponding heading error are taken from the work of Afzal et al. (2010), where indoor magnetic perturbations are studied.It was found that the obtained disturbance on the heading estimate looks similar to a Gaussian curve, and in this analysis it is thus modeled as such.
The disturbance on the heading estimate in time d(t) is modeled as: Here, the amplitude of the disturbance (in radians) is given by A d , the parameter ε controls the width of the Gaussian curve, and t 0 controls the location of the curve in time.For this study, ε = 1, resulting in a disturbance of approximately 4 seconds, and t 0 = 5 s, such that the disturbance occurs at around 5 seconds into the flight.How such a disturbance looks is presented in Fig. 13 for an amplitude A d of 1 rad.
Several amplitudes of the disturbance are tested, namely 0, 0.25, 0.5, 1, and 1.5 rad.The final amplitude of 1.5 results in a maximum heading estimate error of almost 85 • , which is approximately equal to the amplitude of the disturbance shown by Afzal et al. (2010).Note that the disturbance is introduced directly on the measurement of ∆ ψ (the difference in headings between two MAVs).This is the situation that would occur if one of two would fly in a locally perturbed area.
Since the parameter of interest is how the filter for ∑ B compares to the filter for ∑ A , the results are represented as a percentage comparison of the relative localization errors between the two filters.This is visually presented in Fig. 14.In the figure, a positive % means that the filter for ∑ B performs worse than the filter for ∑ A .At 0%, marked by a dotted line, both filters perform equally well.
The comparison shows that as the applied disturbance amplitude on the heading measurement provided to system ∑ A is increased, the region for which ∑ B performs better than ∑ A expands.In the case of the largest disturbance, with A d equal to 1.5 radians, filter ∑ B even performs better at a range noise σ R equal to 8.This result reinforces the presumption that it is not always better to include a heading measurement in the filter, provided that the range measurement is of a high enough accuracy.We will use this insight for the real-world implementation.In the experimental set-up in Sect.4, we will use Ultra Wide Band (UWB) radio modules to obtain range measurements between MAVs.To give an idea of what type of range noise standard deviations can actually be achieved in practice, in the executed experiments with real MAVs, the UWB modules resulted in ranging errors with standard deviations between 0.1 and 0.3.If we assume a normally distributed ranging error, based on the results hown in Fig. 14, it is then clear that the heading-independent system ∑ B would be the preferred choice for all heading disturbance amplitudes (except, trivially, for the situation where there is little to no heading disturbance at all).
4 Leader-follower flight experiment In this section we demonstrate the heading-independent filter in practice, which is used for leader-follower flight in an indoor scenario.

Leader-follower flight considerations
Before designing an actual control method to accomplish leader-follower flight, let's first reflect on the previous observability analysis results from Sect. 2 and their implications with respect to leader-follower flight.We know that in order to have an observable, heading-independent, system, the combined motion of the leader and follower has to meet the observability condition presented in Eq. 33.We further know that in order to to meet this condition, the three intuitive conditions presented by Eq. 34 to Eq. 36 certainly have to be met.Let's first consider these conditions: 1.The first condition (Eq.34) specifies that the relative position between leader and follower must be non-zero.This condition has little implication to leader-follower flight, other than the fact that the follower must follow the leader at a non-zero horizontal distance, which typically is the objective.
2. The second conditions (Eq.35) tells us that both MAVs must be moving.As far as leader-follower flight is concerned, this is automatically accomplished as long as the leader is not stationary.3. The third condition (Eq.36) is especially impactful for leader-follower flight.It specifies that the MAVs should not be moving in parallel (regardless of speed), unless they are also accelerating.A lot of research on leaderfollower flight aims to design control laws that would result in fixed geometrical formations between different agents in the formation.This is typically achieved by specifying desired formation shapes, or desired interagent distances for members in the swarm (Turpin et al., 2012;Gu et al., 2006;Chiew et al., 2015;Saska et al., 2014).By the very nature of fixed geometries, that would result in parallel velocity vectors.
The third condition requires a different approach to leaderfollower flight.Rather than flying in a fixed formation, it is also possible for the follower to fly a delayed version of the leader's trajectory.As long as the leader's trajectory is not a pure straight line for long periods of time, this will result in relative motion between the leader and follower.This is the approach taken in this paper.This solution should also help to prevent the MAVs from getting stuck in an unobservable situation that is not covered by Eq. 34 to Eq. 36, but that is covered by the full observability condition in Eq. 33.We concluded that for the scenarios that are numerically found to be unobservable according to Eq. 33, changing the relative position p only slightly can already result in an observable situation.In the proposed method of having the follower fly a time-delayed version of the leader's trajectory, the relative position vector p will naturally change if the leader's trajectory is not a straight line.

Leader-follower formation control design
We want to construct a leader-follower control method that results in the follower flying a delayed version of the leader's trajectory.As it turns out, this type of control can be directly accomplished with the information provided by the relative localization filter.
Consider the schematic in Fig. 15.It shows two arbitrary trajectories in dotted lines.At the top, in blue, is the trajectory for MAV 1, which is represented by its position vector in time p 1 (t).On the bottom, in orange, is the trajectory for MAV 2, p 2 (t).Suppose the desire is for the follower (MAV 1) to follow the leader's trajectory (MAV 2) with a time delay τ.The control problem for MAV 1 can be expressed as the desire to accomplish p 1 (t) = p 2 (t − τ).
Let t n indicate the current time at which a control input must be calculated.At the current time, MAV 1 has a body Fig. 15 Control problem for leader-follower flight.In blue is MAV 1's trajectory in time p 1 (t).In orange is MAV 2's trajectory in time p 2 (t).The desire is for MAV 1 to drive e(t) to 0 for t → ∞.
fixed reference frame H 1 (t n ), whose origin is p 1 (t n ).At time t n − τ, MAV 1 knows the relative position of the leader in its own body fixed frame H 1 (t n − τ), since this information is provided by the relative localization filter.However, for this control method to work, MAV 1 must have knowledge of where the leader's old position is at the current time t n .This value of interest is depicted by the vector e(t n ) in Fig. 15; it is the positional error with respect to the desired follower's position at time t n .
Let R H i (t 1 )H i (t 2 ) be the rotation matrix from frame H i at time t 2 , to frame H i at time t 1 , defined as: ∆ ψ i | t 2 t 1 is the change in heading angle for MAV i from time t 1 to time t 2 , which can be calculated as: The current positional error for the follower MAV 1, depicted in Fig. 15, can be defined as: The vector ∆ p t n t n −τ represents how much the follower has moved from time t n − τ until t n as defined in frame H 1 (t n − τ).This vector can be calculated using information available to the follower: Finally, one more piece of information is needed in order to be able to design a control law for the follower MAV, which is the model of the follower MAV and how it responds to control inputs.In this paper, it is assumed that the MAV already has stable inner loop control running on board, such that the MAV becomes an outer loop control system that directly can take velocity commands.It is further assumed that with the inner loops in place, the MAV responds like a very simple first order delay filter to velocity commands, such that the differential equation for the its velocity becomes: Where τ τ τ −1 is a diagonal matrix with on the diagonal the inverse values of the time constants that characterize the delay of the system with respect to a control input v 1c .This is only an approximation of how the actual MAV behaves, but it will be shown to be sufficient to accomplish the desired behavior.
With all this information in place, a control law can be designed.The control law is designed using Nonlinear Dynamic Inversion (NDI) principles.In order to use NDI, a state space model is required for the situation at hand.A very similar state space model to the one used for the relative localization filter can be used.Define the state vector as: The state vector is similar to the one defined before for the relative localization filter, with a few small changes.First of all, e = e(t) represents the current positional error for the follower MAV 1 with respect to the leader's old position.Secondly, ∆ ψ and v2 represent again the difference in heading between two MAVs and the velocity of MAV 2, except now ∆ ψ is the difference in heading between frame H 1 (t) and H 2 (t −τ), and v2 is the delayed leader's velocity at time t − τ, such that v2 = v 2 (t − τ).
Similarly, define a new input vector as: Where v 1c is the actual control input fed to MAV 1, and ā2 and r2 represent the same values as a 2 and r 2 , except delayed versions thereof.Therefore ā2 = a 2 (t −τ) and r2 = r 2 (t −τ).Finally, a new set of state differential equations can be defined as: Where R = R(∆ ψ) and S2 = S 2 ( r2 ).
The state that we wish to control is the current positional error that MAV 1 has with respect to the delayed leader's position, so the state e.This state can be represented as: With H given by: The derivative of the control variable with respect to time is equal to: The second derivative of the control variable: With D equal to: This can further be reduced to: At this point the following control law can be chosen: with i now a virtual control input.This control law results in a fully linearized differential equation for the positional error of the follower, since substitution of the control law from Eq. 55 in Eq. 51 results in the following differential equation: Which can be shown to be exponentially stable if the following virtual control is implemented: One of the main findings in the observability study and the simulation results is that the localization error scales more steeply with range noise for system ∑ B than for ∑ A .It is therefore important to use sensors that can provide accurate relative ranging measurements.
In this work, we chose to use Ultra Wide Band (UWB) based radio transceivers.UWB has recently gained attention within the domain of ranging.UWB signals are characterized by their fine temporal and spatial resolution (Correal et al., 2003), which leads UWB based systems to be able to, for example, resolve multipath effects more easily (Win and Scholtz, 1998).Ultimately, this leads to an accurate ranging performance, which is important if using the heading independent filter.Another advantage of UWB is its relative robustness to interference from other radio technologies due to the fact that it operates on an (ultra) wide range of frequencies (Liu et al., 2007;Foerster et al., 2001;Molisch et al., 2006).
The UWB ranging hardware used in the experiments is the ScenSor DWM1000 module sold by Decawave. 2The ranging algorithm that is employed is a particular implementation of the Two-Way Ranging (TWR) method (Neirynck et al., 2016).In order to fuse ranging data with velocity, acceleration, height, and yaw rate data in the localization filter, these variables are also communicated between MAVs by using the UWB devices.The same UWB messages used in the TWR protocol are also used to communicate these variables.
The UWB module transceiver has been installed on the Parrot Bebop 2 platform. 3The Bebop 2 runs custom autopilot software designed using the open-source autopilot framework Paparazzi UAV. 4 Paparazzi UAV provides the stable inner loop control loops for the Bebop 2 using Incremental NDI (INDI).This allows us to control the outer loop by giving the computed velocity commands to the INDI inner loops.
Velocity and height measurements are also necessary for the relative localization filter.In the initial experiments, they are provided by an overhead Motion Capture System (MCS) by OptiTrack. 5In a second iteration of the experiment, they are fully provided by on-board sensors.The velocity data is obtained from the MAVs' on-board bottom-facing camera using Lucas-Kanade optical flow.Height is measured using an on-board ultrasonic sensor that the Bebop 2 is equipped with by default.At all times, the acceleration and yaw rate measurements are obtained from the MAVs' on-board accelerometers and gyroscope, respectively.The experiments are first conducted with two MAVs (one leader and one follower), detailed in Sect.4.4, and then performed again with three MAVs (one leader and two followers), detailed in Sect.4.5.

Leader-follower flight with one follower
The experiment with one follower MAV consists of one Bebop 2 following another Bebop 2 using the control law presented in Sect.4.2.At first, right after take off, the MAVs fly concentric circles just like the ones shown in Fig. 11.This procedure is there to make sure that the EKF running onboard the MAVs has time to converge to the correct result, such that by the time the follower MAV is instructed to start following the leader, the follower has a correct estimate of the relative location of the leader.
When leader-follower flight is engaged, the trajectory of the leader has been designed to sufficiently excite the the relative localization filter during the leader-follower flight and to decrease the likelihood of being stuck in unobservable states.This has been done by introducing frequent turns in the trajectory to have changing relative velocities and accelerations.The follower is instructed to follow the leader's trajectory with a time delay of τ = 5 seconds.
It is important to note that, for safety reasons, the norm of the follower's commanded velocity ||v 1c || 2 during both experiments is saturated at 1.5 m/s.The measure is taken because the MAVs were flying in a relatively small confined area (10 m by 10 m).This change does however have consequences for the performance of the follower's tracking, which is discussed further in the next sections.

Leader-follower flight with velocity and height information from a MCS
First, the case where velocity and height information is provided by the MCS is studied.In Fig. 16 the trajectory flown by the follower is compared to the trajectory of the leader.The x and y coordinates are compared separately for part of the flight in Fig. 17a and Fig. 17b.In Fig. 18, a time composition of overhead camera images is given for 5 seconds of flight as an illustration.The follower's position is shown at seven time instances during these 5 seconds, and is compared to the leader's trajectory.
A total of 200 seconds of leader-follower flight were logged and will be analyzed here.During this time, several laps of the designed trajectory were executed.The trajectories in Fig. 16 to Fig. 18 indeed show that the follower is successfully tracking a delayed version of the leader's trajectory.The actual error distribution for the norm of the relative location estimate ||p|| 2 is shown in Fig. 19.The errors have a mean value of 18.4 cm and a maximum value of 77.5 cm, at maximum inter-MAV distances up to 5 m.
- 4 -3 -2 -1 0 1 2 3 4 5   Since, in this experiment, the velocity and height measurements were provided with high accuracy by the MCS, one would expect the primary source for the localization error to be the ranging error from the UWB modules.How-  error is caused by a relative localization error from the follower's perspective, which will inevitably affect the tracking performance.However, since the relative localization error is considerably lower than the tracking error, there must be more sources to the error.
One source of error is the fact that the follower's response to a velocity command v 1c is modeled as a first order delay.In reality, the MAV has some overshoot with respect to commands, which is not modeled by this first order delay.This model mismatch by itself might not be that harmful to the performance, since the control law would respond with more aggressive velocity commands as a reaction to the MAV not behaving as modeled.However, the control law's freedom is severely restricted by the command saturation at 1.5 m/s, which means that the follower cannot move as fast as the command law demands.This argument is further supported by a qualitative analysis of the follower's trajectory with respect to the leader's trajectory in Fig. 16.The trajectory of the follower often seems to take 'shortcuts' with respect to the leader's trajectory.This falls in line with the expected behavior due to the command saturation.The control law is designed not only to track the trajectory of the leader in space, but also in time.As the follower starts lagging behind the leader more than the desired τ = 5 seconds, the follower starts to take shortcuts in the trajectory to catch up with the leader.This error would be less prevalent if the command saturation were increased.

Leader-follower flight with only on-board measurements
We now demonstrate the workings of the proposed methods in this paper when only on-board sensing is used.In this setup, the follower MAV does not use any MCS information.Instead, the velocity information comes from Lucas-Kanade optical flow measurements while the height is derived from the on-board ultrasonic sensor.Similarly, the leader MAV directly communicates optical flow velocities and ultrasonic height measurements (along with accelerations and yaw rate -4 -3 -2 -1 0 1 2 3 4 5  Fig. 23 The trajectory of the follower compared to the delayed trajectory of the leader for the experiment with only on-board sensing.from the IMU) to the follower MAV for use in the relative localization filter.The MCS is only used to log ground truth data and for the leader to safely fly its trajectory.No MCS data is used by the follower at all.Again, 200 s of leaderfollower flight with full on-board sensing took place successfully and will be analyzed here.
The trajectory of the follower with respect to the delayed leader's trajectory is compared in Fig. 22 and Fig. 23.Furthermore, another time composition for 5 seconds of flight where the follower is tracking the leader is given in Fig. 24.
The main qualitative difference with respect to the situation where the MCS was still used for velocity and height information is the fact that the follower's trajectory appears less smooth.Otherwise, the performance seems qualitatively similar.The follower still appears to take 'shortcuts' with respect to the leader's trajectory, although the increased disorder in the follower's trajectory makes this less apparent.
The tracking error distribution for the on-board sensing case is given in Fig. 26.The mean tracking error is 50.8 cm and the maximum error is 1.47 m.The relative localization error is given in Fig. 25.Here, the mean error is 22.6 cm and the maximum error is 75.8 cm, at maximum MAV distances up to 5.2 meters.
The performance when using only on-board sensing is very similar to when using the MCS for height and velocity data.This can be mainly attributed to the fact that the measurements that have been replaced (the height and ve-  The primary reason as to why the trajectory of the follower with on-board sensors still seems slightly more disordered is the fact that the follower has difficulty to accurately control its altitude when using only on-board sensing.The follower now purely relies on height measurements from its ultrasonic sensor.The update rate of this sensor is low, and in between measurements the follower uses (noisy) accelerometer data to update its height.This sometimes causes the follower to believe its altitude is different than it really is, causing it to rapidly ascend or descend.This takes up thrust, restricting the follower's ability to maneuverer accurately in the horizontal plane due to thrust saturation.To demonstrate that the methods in this paper can also scale to more than one follower, the leader-follower flight is also performed with two follower MAVs instead of one.This is done both with MCS height and velocity data and with only on-board sensing.
For this purpose, The UWB messaging protocol is adapted to allow every MAV to perform ranging with every other MAV.The MAVs also communicate a unique (pre-assigned) identification number within the UWB messages.The followers can use this identification number to determine which messages originate from the leader so that they individually keep track of the leader as before.The main consequence of the increased messages is a drop in the UWB range update rate, which is reduced from about 25 Hz with 2 MAVs, to about 16 Hz with 3 MAVs.
This time, due to the lack of space available, there is no initialization flight procedure to give the EKFs of the followers time to converge.Instead, the MAVs are placed in starting positions and orientations that roughly match with what EKFs on-board the MAVs are initialized to.Although this placement is done purely by eye, it is proven to be sufficient to safely start the leader-follower flight.
The leader flies the same trajectory as before.The first follower follows this trajectory with a τ = 4 second delay, and the second follower follows it with an τ = 8 second delay.Once again, 200 seconds of successful flight data is logged and analyzed.
An overhead camera image for the flight with MCS height and velocity data is presented in Fig. 27, giving an idea of how the experiment looked like.The trajectories for this flight are displayed in Fig. 28 for the leader and two followers.For the flights with only on-board information, the trajectories are shown in Fig. 29.
As for the case with just one follower, we see that the followers tend to take shortcuts with respect to the leader's trajectory.Furthermore, the flights using only on-board in- Fig. 28 Trajectory of leader and two followers using MCS height and velocity formation are less smooth than those with MCS height and velocity information.For the flight with MCS data, follower 1 has a MAE for the relative localization error of only 15.8 cm.By comparison, follower 2 has a MAE of 43.9 cm.Furthermore, followers 1 and 2 have MAE for the tracking of 42.9 cm and 70.3 cm, respectively.The flight with only onboard sensing resulted in a relative localization MAE of 51.8 and 53.6 cm.The tracking MAE this time was 58.6 and 98.4 cm.

Comparison of flights
In this section we present the relative localization and tracking MAE of the various flights that were executed.We also discuss in more detail the most noteworthy differences between experiments.
All the errors are presented in Tab. 2. The first noteworthy observation is the fact that, for the experiment with two followers, the tracking performance of the second follower is worse than for the first follower in both the MCS and fully on-board case.This is a byproduct of the fact that the pro-  posed leader-follower control method inherently relies on integration of velocity information in time.As the delay with which the follower must follow the leader increases, so does the period of time over which the follower must integrate its velocity.This is subject to drift, which shows in the tracking performance.This effect is more noticeable in the fully onboard case, since the velocity estimates from optical flow methods are less accurate than the ones computed by the MCS.
Another result is that the localization error for follower 2 in the MCS case is higher than for the first follower.This can be explained, in part, by the fact that follower 2 has a larger mean range with respect to the leader than follower 1 does (4.2 m compared to 2.9 m).To inspect this deeper, we looked at the logged range between the MAVs.It was found that follower 2 had substantially larger ranging errors with the leader than follower 1.This can be appreciated in Fig. 30, where the ranging error distributions are compared.In both cases, the mean is close to zero, yet the distribution for follower 2 is significantly wider.At this stage, it is not clear what the primary cause for this drop in range error is.It shall be studied further in future work when the scalability of the system is addressed in more detail.
A final result that stands out is that both followers 1 and 2 have substantially higher localization errors in the onboard case than was found for the on-board experiment with a single follower.This result appears to be due to a combination of factors.The increased communication traffic caused a decrease in the filter update rate and also resulted in an increase in ranging frames dropped.Follower 2, as mentioned above, showed a worse ranging performance than follower 1. Follower 1, in turn, had slightly less accurate optical flow velocity estimates than were obtained with the single follower flight (21 cm/s MAE compared to 15 cm/s before) and also slightly higher ranging errors than for the single follower flight (15 cm MAE compared to 8 cm before).All factors combined, both followers suffered a comparable degradation in localization performance.

Discussion
In this section we revisit the observability analysis from Sect. 2 with the obtained experimental data.We also present some remarks on the scalability of this methodology to larger groups of MAVs.

Remarks on observability
In Sect.2.5 showed that for a specific set of velocity, accelerations and relative positions for both MAVs, the system will become unobservable.To directly integrate the full observability condition in the design of a leader-follower system is difficult due to its high dimensionality.By having followers fly a delayed version of the leader's trajectory, it is possible to naturally vary the relative positions between leader and follower, as long as the leader's velocity changes in time.
Given the sparsity of unobservable relative positions, we therefore postulated that this control behavior would be sufficient to limit unobservable situations.Furthermore, even if an unobservable situation were to occur, this would only be for a short period of time, as the relative position continuously changes and the system automatically transitions back to being observable.
Having performed the experiments and collected all the ground truth data, it is now possible to test whether this assumption is valid.All the parameters needed to evaluate Eq. 33 have been logged during the experiments and can be inserted into Eq.33 to check the observability of the relative localization filter in time.In line with our previous analysis, the measure of observability of the system is represented by the cross product between the left hand side of Eq. 33 and the relative position vector p.Once more, we shall take a threshold of 1, meaning that an observability value between -1 and 1 is considered unobservable.Although theoretically only a value of 0 would indicate an unobservable system, the higher threshold is chosen to account for noise in the data.
With the chosen threshold, the unobservable data points for the MCS and the on-board flight are 4.76% and 4.75% of all the data points, respectively.The unobservable points are spread in time, thus giving the system ample observable data in between to recover from the short periods of unobservability.Furthermore, isolated events of unobservability are not expected to cause issues.Instead, they can gradually cause an increase in the localization error in time.This has also been confirmed by the simulations in Sect.3.
Further qualitative inspection of the data does not show a correlation between the unobservable regions of the flight and the relative localization error.To demonstrate this, the localization error is compared to the observability of the filter in Fig. 31 for a small segment of the flight with MCS information.For easier comparison, the observability has been reduced to a binary value, where a value of '1' indicates that the system is within the threshold of unobservability at that time.It can be seen that there is no apparent correlation between the two parameters.

Remarks on scalability
The experimental results in Sect. 4 show that the methods in this paper can successfully scale to two followers that follow a leader in a confined area.Even when full on-board sensing is used by the followers, more than three minutes of successful autonomous flight were demonstrated, with no pilot input.
Despite the successful results, analysis of the data does show a substantial rise in localization and tracking errors when scaling up to two MAVs.This raises the question of what would happen if even more MAVs are added to the experiment; would this be viable?
One of the results we found is that there is a correlation between the tracking performance of the follower and the time delay with which it follows the leader's trajectory.The follower that tracked with a time-delay of eight seconds showed consistently larger tracking errors than the followers with four or five second delays.An alternative solution to the two follower problem is having one follower follow the leader and the other following the first follower.With such an arrangement, both followers could follow another MAV with the same time delay.This setup has not yet been studied in this work, but could prove to be a better alternative to explore in future research.
Another result we found is that the update rate reduces when flying with two followers instead of one.It is to be expected that adding more MAVs requires additional data communication, yet a drop from 25 Hz to 16 Hz is quite significant for adding just one more MAV.The main remark to make here is that this reduction in update rate is very much dependent on the software and hardware used for these experiments.It should be possible to significantly increase the update rate to allow for more MAVs without sacrificing the update rate to a large extent.
As an example, in these experiments we operated the UWB modules on the lowest data rate settings (110 kbps).Furthermore, every message contains a lengthy preamble of 2048 bits, resulting in substantial protocol overhead for every transmitted message (the actual payload of the UWB messages is less than 200 bits).This should theoretically help to improve the ranging accuracy, but in practice will most likely not make a big difference at the small inter-MAV ranges occurring in these experiments (DecaWave, 2017).The maximum data rate that the UWB modules support is actually 6.8 Mbps and the preamble can be as short as 64 bits.These would allow for much higher update rates, even with three or more MAVs.One would, however, need to examine what such a change would have on ranging accuracy and stability.

Conclusion
The work in this paper has shown the feasibility of headingindependent range-based relative localization on MAVs.We now know that removing the dependency on a common heading between MAVs has two main disadvantages: the motion of agents must meet more stringent conditions to be observable and the relative localization becomes more susceptible to noise on the range measurements.The clear advantage, on the other hand, is that the filter is no longer affected by local disturbances in Earth's magnetic field.As shown by our simulations, small magnetic perturbations can already lead to a large negative impact, showing how a headingindependent method can actually perform better than the heading-dependent method.
The results of our observability analysis have shown that leader-follower flight is a difficult task when using the proposed relative localization method.Fixed geometry formation flight is not possible.Instead, we developed a method that allows one MAV to follow another MAV's trajectory with a certain time delay.This approach has been shown to stay sufficiently clear from unobservable conditions, which has allowed us to successfully demonstrate leader-follower flight in practice.
Using only on-board sensory information, one MAV can localize another MAV with a mean error of just 22.6 cm over 200 seconds of leader-follower flight.This consequently allows the MAV to track another MAV's trajectory with a mean error of 50.8 cm.The method has been demonstrated to work also with two followers tracking the same leader.

Future work
There are plenty of opportunities to research within the domain of range based relative localization.Certainly, one such opportunity is the initial convergence behavior of the filter.The initial estimate of the EKF is important to quickly converge to a correct estimate of the relative location of another MAV.If the initial condition is too different from the real situation, the filter have difficulties to converge.One primary problem is that there exist ambiguous states where the EKF can converge to and from which it is difficult to then escape.In the future, it would thus be interesting to research methods to address this problem.Examples solutions could be alternative filters (e.g., a particle filter), or running multiple filters in known ambiguous states to identify the correct state more easily.
Furthermore, the current leader-follower implementation flight uses a large amount of past data values, and directly uses state values like the velocities of the two MAVs to implement its control method.It would be interesting to research other methods of accomplishing this type of leaderfollower flight.For example, it might be possible to perform real time polynomial data fitting on the relative positions of the leader.The resulting polynomial trajectories could instead be used to obtain the velocities and accelerations through analytical derivations of the polynomials.This might result in less data that needs to be stored and smoother trajectories.
Finally, considering the hardware used in the experiments, the importance of consistent, high frequency communication and ranging has become apparent.It would be valu-

Fig. 4
Fig. 4 Color map of observability for different relative positions.The velocities and accelerations of the MAVs are kept as depicted by figure 3 and the values for p = [p x , p y ] are varied over a 10 m range.

Fig. 13
Fig. 13 Disturbance on the relative heading measurement in time, for an amplitude A d of 1 rad

Fig. 14
Fig. 14 Percentage error comparison between ∑ B and ∑ A for different disturbance amplitudes A d .Positive percentage means ∑ B performs worse than ∑ A .

Fig. 16
Fig.16The trajectories of leader and follower during experiment with MCS height and velocity

Fig. 18
Fig. 18 Time composition of overhead camera images of leader and follower MAV in time, for the experiment with MCS height and velocity.Indicated in orange and marked by p 2 (t), is part of the leader's trajectory.The leader's final position is indicated by p 2 (t = 0).Seven points in time of the follower's trajectory are indicated in the image.According to the control objective, p 1 (t = 5) should equal p 2 (t = 0).

Fig. 19
Fig. 19 Histogram of the localization error for the follower during experiment with MCS height and velocity

Fig. 21
Fig. 21 Histogram of the tracking error ||e|| 2 for the follower during experiment with MCS height and velocity

Fig. 22
Fig. 22 Trajectory of leader and follower during experiment with only on-board sensing and processing

Fig. 24
Fig. 24 Time composition of overhead camera images of leader and follower MAV in time, for the experiment with only on-board sensing.Indicated in orange and marked by p 2 (t), is part of the leader's trajectory.The leader's final position is indicated by p 2 (t = 0).Six points in time of the follower's trajectory are indicated in the image.According to the control objective, p 1 (t = 5) should equal p 2 (t = 0).

Fig. 25
Fig. 25 Histogram of the localization error for the follower during experiment with only on-board sensing and processing

Fig. 26
Fig. 26 Histogram of the tracking error ||e|| 2 for the follower during experiment with only on-board sensing and processing

Fig. 27
Fig.27Overhead camera image of leader and two followers using MCS height and velocity.In orange is the leader's trajectory marked at 0.5 second intervals.

Fig. 29
Fig. 29 Trajectory of leader and two followers using only on-board information

Fig. 30
Fig. 30 Comparison between ranging error distributions for follower 1 and 2 for the flight with MCS height and velocity data.

Fig. 31
Fig.31Comparison between localization error and the observability of the filter.An unobservable value of '1' means the observability measure is within the threshold of unobservability (between -1 and 1)

Table 1
Average Mean Absolute Error for ∑ A and ∑ B over 1000 runs with different noise standard deviation on the range measurement 2 Fig. 12 AMAE in estimate of ||p|| 2 for ∑ B and ∑ A

Table 2
Comparison of mean localization (loc.)errors and mean tracking (track.)errors for all performed experimental flights, both for MCS and fully on-board (on-b.)flights.