Moses Bangura
email: mbangura@andrew.cmu.edu

Xiaolei Hou
email: hou.xiaolei@nwpu.edu.cn

Guillaume Allibert
email: allibert@i3s.unice.fr

Robert Mahony
email: robert.mahony@anu.edu.au

Nathan Michael
email: nmichael@andrew.cmu.edu

Supervisory Control of Multirotor Vehicles in Challenging Conditions using Inertial Measurements

We consider the problem where a supervisor or remote pilot provides a real-time linear velocity reference to a multirotor aerial robot; either through a traditional remote control handset, a modern haptic interface, or semi-autonomous guidance control system. In all such cases, the goal is to servocontrol the vehicle's velocity to the set point as quickly and as efficiently as possible. The challenge is to achieve this robustly in the presence of unknown wind disturbances and in situations where the vehicle moves into GPS denied environments (indoors, urban canyons, forests) where estimation of the vehicle's velocity is challenging. These situations include unclutterred environments, poor visibility environments caused by poor lighting and poorly textured visual environments where laser and vision based sensors become unreliable. The approach taken is to develop a coupled non-linear complementary velocity aided attitude filter that provides estimates of both the inertial and body-fixed frame linear velocities, as well as the attitude of a multirotor aerial vehicle, that functions effectively even when only the inertial measurement unit (IMU) and barometric sensor measurements are available. When full inertial velocity measurements are available (from GPS, Vicon or a vision system), the filter additionally estimates the external wind speed. In this paper we formally present the proposed filter along with experimental results and a comparison of the filter to recent results in the literature and in situations where inertial reference frame velocities are available intermittently. The proposed filter is computationally simple to implement and easy to calibrate, tune and provides an excellent base level functionality for modern multirotor aerial robotic systems that will be required to function robustly in a variety of environments.

I. INTRODUCTION

Multirotors are aerial vehicles with multiple pairs of motorrotor systems for the generation of thrust and vehicle dynamic control [START_REF] Bangura | Nonlinear dynamic modeling for high performance control of a quadrotor[END_REF]. Due to their light weight, simple dynamics, cost and availability, they have become the standard platform for both aerial robotics research and the hobby drone communities. Multirotor robots have numerous potential military and civilian applications including construction of structures [START_REF] Lindsey | Construction of cubic structures with quadrotor teams[END_REF], search and rescue [START_REF] Waharte | Supporting search and rescue operations with uavs[END_REF], surveillance of earthquake buildings [START_REF] Michael | Collaborative mapping of an earthquake-damaged building via ground and aerial robots[END_REF] and many more. Some of these environments are indoor or occluded from global position systems (GPS) and require a secondary sensor modality to estimate linear velocity. Vision systems are a natural solution to this challenge and are presently a hot topic in the robotics and computer vision community. Vision based approaches have led to visual inertial navigation system (VINS) based state estimation gaining increasing popularity in the community [START_REF] Shen | Vision-based state estimation and trajectory control towards high-speed flight with a quadrotor[END_REF], [START_REF] Faessler | Autonomous, vision-based flight and live dense 3d mapping with a quadrotor micro aerial vehicle[END_REF], [START_REF] Mebarki | Nonlinear visual control of unmanned aerial vehicles in gps-denied environments[END_REF]. These approaches generally use an unscented Kalman filter (UKF) [START_REF] Shen | Vision-based state estimation and trajectory control towards high-speed flight with a quadrotor[END_REF], extended Kalman filter (EKF), Gaussian particle filter [START_REF] Bry | Aggressive flight of fixed-wing and quadrotor aircraft in dense indoor environments[END_REF] or a sliding window observer [START_REF] Shen | Initializationfree monocular visual-inertial state estimation with application to autonomous mavs[END_REF] based on a combination of lasers, lidars, RGB and RGB-D cameras, and stereovision. An example of these algorithms that has been used to do indoor teleoperation is direct visual odometry [START_REF] Stegagno | A semiautonomous uav platform for indoor remote operation with visual and haptic feedback[END_REF] using an RGB-D sensor. Table I provides a summary of some of the current vision, laser and LiDAR based techniques being used for localisation in GPS denied environments. These techniques however, are computationally challenging and aerial robots with these sensing modalities often carry secondary CPU boards such as Intel Atom board [START_REF] Grzonka | A fully autonomous indoor quadrotor[END_REF], Odroid XU4 or a dedicated Intel Core i7 ground station [START_REF] Faessler | Autonomous, vision-based flight and live dense 3d mapping with a quadrotor micro aerial vehicle[END_REF]. Even the most computationally tractable vision based algorithm, optical flow, which estimates the ego-motion of the vehicle as scaled linear and angular velocities require significant computational resources [START_REF] Grabe | On-board velocity estimation and closed-loop control of a quadrotor uav based on optical flow[END_REF], [START_REF] Grabe | Robust optical-flow based self-motion estimation for a quadrotor uav[END_REF], [START_REF] Grabe | Nonlinear ego-motion estimation from optical flow for online control of a quadrotor uav[END_REF], [START_REF] Briod | Optic-flow based control of a 46g quadrotor[END_REF]. It is expected then that most aerial robotic systems will not carry multiple vision systems, and should their single vision system fail or be obscured by fog, rain, dust or low visibility due to low light or poor texture, or should the environment they are observing be highly dynamic, invalidating the static environment assumption that is fundamental in most VINS, then the functionality of the system would be significantly compromised. This is particularly important in the case of supervisory control activities where a human is integral in task and such situations may well be encountered. For supervisory control linear velocity reference, derived from user inputs are provided in real-time to the vehicle [START_REF] Lam | Artificial force field for haptic feedback in uav teleoperation[END_REF], [START_REF] Mersha | On bilateral teleoperation of aerial robots[END_REF], [START_REF] Hou | Comparative study of haptic interfaces for bilateral teleoperation of vtol aerial robots[END_REF]. Although some supervisory systems are based on trajectories [START_REF] Yang | A framework for efficient teleoperation via online adaptation[END_REF] and position mapping [START_REF] Farkhatdinov | Hybrid position-position and positionspeed command strategy for the bilateral teleoperation of a mobile robot[END_REF], in all cases a good onboard estimate of the vehicle velocity is critical to implementation of a robust and reliable supervisory control. Thus, it is critical that an aerial robot has robust estimates of its linear velocity available at all times, including when GPS is unavailable or vision systems fail.

Direct airspeed measurements using pitot tubes [START_REF] Yeo | Onboard flow sensing for downwash detection and avoidance with a small quadrotor helicopter[END_REF], [START_REF] Arain | Real-time wind speed estimation and compensation for improved flight[END_REF] has been considered for quadrotor vehicles. This method suffers from accuracy and slow dynamic response of pitot tube technology; typically around 100ms rise time. Arain et. al. [START_REF] Arain | Real-time wind speed estimation and compensation for improved flight[END_REF] noted that pitot tube velocity measurements for quadrotors is unreliable for velocities under 1m/s and they observed errors of up to 2m/s for ground truth forward velocity of 6m/s. Yeo et. al. [START_REF] Yeo | Onboard flow sensing for downwash detection and avoidance with a small quadrotor helicopter[END_REF] used four pitot tubes mounted underneath

Class of Algorithm Measurements Comments Optical Flow/Visual Odometry ∆x, ∆ȳ, z

These algorithms require a camera and can be used along with an IMU to obtain de-rotated optical flow measurements. These measurements can be scaled by range sensors to obtain twists (position change between successive measurements) ∆ζ in x and y. Example algorithms are the PX4FLOW [START_REF] Grabe | On-board velocity estimation and closed-loop control of a quadrotor uav based on optical flow[END_REF]. Laser odometry ∆x, ∆ȳ, z LiDARs can be used to estimate the twist in the axis of measurement. Example algorithm is the RF2O laser odometry algorithm [START_REF] Jaimez | Planar odometry from a radial laser scanner. a range flow-based approach[END_REF].

Laser SLAM

x, ȳ, z Uses laser scan matching of LiDAR measurements. The algorithms are usually accurate to generally not require loop closures. Example algorithm is Hector mapping [START_REF] Kohlbrecher | A flexible and scalable slam system with full 3d motion estimation[END_REF].

Monocular

Visual SLAM

x, ȳ, z

Uses a camera to determine the scaled position and exact rotation of a camera using SIFT/SURF features. In addition, these algorithms require an IMU and additional estimation mechanism to estimate the absolute scale in order to estimate the absolute pose from a camera image [START_REF] Shen | Initializationfree monocular visual-inertial state estimation with application to autonomous mavs[END_REF]. Examples include SVO, ORB-SLAM, LSD.

RGBD

Odometry/SLAM

x, ȳ, z

It uses depth information (x, y, z) of observed points in an RGB image to estimate the position and attitude of the sensor. The range and accuracy of the algorithm depend on the sensor range and noise characteristics [START_REF] Mur-Artal | Orb-slam: a versatile and accurate monocular slam system[END_REF].

Stereo SLAM

x, ȳ, z This is the most computationally intensive SLAM algorithm and uses stereo images to do 6DoF SLAM [START_REF] Mur-Artal | Orb-slam: a versatile and accurate monocular slam system[END_REF].

each rotor to measure the axial velocities through the rotors. They also noted slow rise time and obtained errors of up to 0.4m/s for air speed of 1.5m/s. Tomić et. al. [START_REF] Tomić | The flying anemometer: Unified estimation of wind velocity from aerodynamic power and wrenches[END_REF] used aerodynamic power [START_REF] Bangura | Aerodynamic power control for multirotor aerial vehicles[END_REF] to estimate induced flow through the rotor, however, there is insufficient information with this approach to discriminate between motions in e 1 , e 2 and e 3 directions [START_REF] Bangura | Thrust control for multirotor aerial vehicles[END_REF]. Davis and Pounds [START_REF] Davis | Direct sensing of thrust and velocity for a quadrotor rotor array[END_REF] used a force torque sensor to estimate total velocities of up to 1m/s with an accuracy of about 0.1m/s. A more promising approach is based on exploiting the linear induced drag effects of rotor motion [START_REF] Martin | The true role of accelerometer feedback in quadrotor control[END_REF] to estimate translational linear velocities. This drag model along with a barometer was used in an extended Kalman filter to estimate all three linear velocities [START_REF] Burri | Robust state estimation for micro aerial vehicles based on system dynamics[END_REF]. Furthermore, Sikkel et. al. [START_REF] Sikkel | A novel online model-based wind estimation approach for quadrotor micro air vehicles using low cost mems imus[END_REF] used this model along with GPS to develop a nonlinear extended Kalman filter for predicting wind speed.

State-estimators for quadrotors can be grouped into attitude, velocity and position and the more recent combined estimation for attitude and linear velocity, velocity aided attitude observers. The most common of these is the attitude observer which include the left-invariant extended Kalman filter (LIEKF) [START_REF] Bonnabel | Left-invariant extended kalman filter and attitude estimation[END_REF], [START_REF] Barrau | The invariant extended kalman filter as a stable observer[END_REF], and non-linear complementary filters [START_REF] Mahony | Nonlinear complementary filters on the special orthogonal group[END_REF]. The major disadvantage of attitude only estimation is the requirement for a second velocity and/or position filter. Velocity aided attitude estimation combines these two state estimation problems into a single observation problem. The vehicle's linear velocity measurements can be posed in either the inertial frame [START_REF] Hua | Attitude estimation for accelerated vehicles using gps/ins measurements[END_REF], [START_REF] Roberts | On the attitude estimation of accelerating rigid-bodies using gps and imu measurements[END_REF], [START_REF] Martin | An invariant observer for earth-velocity-aided attitude heading reference systems[END_REF] or in the body fixed frame [START_REF] Allibert | Estimating body-fixed frame velocity and attitude from inertial measurements for a quadrotor vehicle[END_REF], [START_REF] Allibert | Velocity aided attitude estimation for aerial robotic vehicles using latent rotation scaling[END_REF], [START_REF] Hua | Stability analysis of velocityaided attitude observers for accelerated vehicles[END_REF]. Bonnabel [START_REF] Bonnabel | Invariant extended kalman filter: theory and application to a velocity-aided attitude estimation problem[END_REF] proposed an invariant extended Kalman filter (EKF) for velocity and attitude estimation in quaternions. Abeywardena et. al. [START_REF] Abeywardena | Improved state estimation in quadrotor mavs: A novel drift-free velocity estimator[END_REF] and Leishman et. al. [START_REF] Leishman | Quadrotors and accelerometers: State estimation with an improved dynamic model[END_REF] used an Euler angle based combined roll and pitch and translational velocities with an extended Kalman filter (EKF) formulation. Hua et. al. developed nonlinear observers based on Lyapunov analysis and approximation arguments [START_REF] Hua | Attitude estimation for accelerated vehicles using gps/ins measurements[END_REF], [START_REF] Hua | Velocity-aided attitude estimation for accelerated rigid bodies[END_REF]. The authors prior work [START_REF] Allibert | Estimating body-fixed frame velocity and attitude from inertial measurements for a quadrotor vehicle[END_REF], [START_REF] Allibert | Velocity aided attitude estimation for aerial robotic vehicles using latent rotation scaling[END_REF] also considers the problem from a nonlinear observer design perspective in the body-fixed frame.

In this paper we make the following contributions:

1) We propose a non-linear state-estimation and control algorithm for supervisory control of aerial vehicles that uses all available sensors systems (GPS, vision, IMU, etc) but is robust to loss of all sensors except the proprioceptive IMU and barometer sensor systems. 2) We provide experimental verification of the proposed closed-loop performance and comparison to state-of-theart algorithms. 3) We demonstrate the robustness of the closed-loop system to loss of exteroceptive sensors systems (GPS/Vision/etc).

This work builds on our previous work in state estimation and control of aerial vehicles [START_REF] Bangura | Non-linear velocity aided attitude estimation and velocity control for quadrotors[END_REF], [START_REF] Allibert | Estimating body-fixed frame velocity and attitude from inertial measurements for a quadrotor vehicle[END_REF], [START_REF] Allibert | Velocity aided attitude estimation for aerial robotic vehicles using latent rotation scaling[END_REF]. A key innovation of the approach is to run a coupled pair of non-linear velocity aided attitude filters, one that estimates the body-fixed frame velocity, and one that estimates the inertial velocity of the vehicle. The body-fixed frame filter uses the inertial measurement unit (IMU), magnetometer and barometer along with a drag model to estimate attitude and linear velocity expressed in the body-fixed frame. This filter requires only proprioceptive sensor systems on board the vehicle and is robust to GPS drop-out, loss of visual texture or failure of vision systems due to loss of visual texture or presence of dust or other visual impediments. However, the filter depends on velocity estimates derived from a dynamic model of rotor drag and based on noisy accelerometer measurements, as well as the barometer sensor that is subject to variable bias. As a consequence, the body-fixed velocity aided attitude filter does not have the same performance as filters that use more precise sensors such as GPS or vision INS. The inertial frame filter uses the IMU in the angular velocity estimation along with inertial reference velocity measurements derived from GPS, vision, or other exteroceptive sensors. This filter has a high performance, providing precise state-estimates, however, it fails should the exteroceptive sensor systems fail. The proposed coupled filter architecture synchronises the two filter state estimates and ensures that the high performance aspects of the inertialframe velocity aided attitude filter are exploited when they are available, and that in the event that the exteroceptive sensor systems fail, the body-frame velocity aided attitude filter will continue to function and keep the inertial frame filter synchronised so it is available to come on line the moment the GPS or vision INS signal is available. An added advantage of the approach is that we obtain an estimate of the average wind conditions in which the vehicle is flying from the constant offset estimate between linear velocity estimates (transformed into the same frame of reference). A teleoperation scheme with a human operator is used to demonstrate the proposed supervisory control algorithm, although the approach is applicable to any real-time velocity reference signal. The proposed control ensures locally exponentially stable velocity tracking.

Experimental results are presented to show performance of the scheme both with and without linear inertial reference frame velocity measurements. Results are also presented to demonstrate the computational efficiency of the algorithm, accuracy and superiority to recent results obtained in the literature.

The remainder of the paper is organised as follows: the nonlinear quadrotor model and the means of obtaining the measurements and estimates of these measurements is described in Section II, the control of the vehicle using the estimated velocities and attitude is described in Section III, in Section IV, we present the hardware system and in Section V, we present experimental results and comparisons to the state-of-the art algorithms.

II. OBTAINING MEASUREMENTS AND STATE ESTIMATION

In this section, we present the non-linear quadrotor model in freestream with no external forces written in both the body-fixed and inertial reference frames. These models form the bases of the state estimation scheme presented later in the section (Section II-D). The model is also used in the the development of the velocity controller in Section III-A. The different sensors that can be used to obtain the various measurements of the different variables are also described. The measurements are then used in our proposed coupled nonlinear complementary velocity aided attitude filter.

In the sequel, e 1 , e 2 , e 3 ∈ R 3 will be used to denote unit vectors in x, y and z directions respectively.

A. Quadrotor modelling

To present the non-linear model, consider Figure 1 which has an inertial-fixed frame denoted by {A} and a body-fixed frame {B}. If the mass of the quadrotor is m, g, acceleration due to gravity and T is the total thrust or heave force, v ∈ R 3 is the linear velocity of the vehicle in {A} and R ∈ SO(3) is the rotation matrix from {B} to {A}, D ∈ R 3 is the drag force expressed in {B}, and Ω ∈ R 3 is the angular velocity of the vehicle in inertial frame, then the linear dynamics model in {A} is given by [

1] v = g e 3 - T m R e 3 + 1 m RD, (1a)
Ṙ = RΩ × , (1b)
where of the vehicle in the body-fixed frame, then (1) expressed in {B} is

Ω × ∈ R 3×3 is the skew-symmetric matrix. It is such that Ω × w = Ω×w for all w ∈ R 3 . If V = R ⊤ v is the velocity
V = -Ω × V + gR ⊤ e 3 - T m e 3 + D m , (2a)
Ṙ = RΩ × . (2b)
In the sequel, if x is the true value of a quantity, then the following notation x, x, A x, B x are used to denote its measurement, estimated value expressed in reference frames A and B respectively.

B. Body-fixed frame measurement of linear velocity B V

We take an approach based on exploiting the linear induced drag force generated by rotors to estimate the horizontal velocity displacement of the vehicle [START_REF] Martin | The true role of accelerometer feedback in quadrotor control[END_REF], [START_REF] Hua | Attitude estimation for accelerated vehicles using gps/ins measurements[END_REF], [START_REF] Abeywardena | Improved state estimation in quadrotor mavs: A novel drift-free velocity estimator[END_REF], [START_REF] Hua | Velocity-aided attitude estimation for accelerated rigid bodies[END_REF], [START_REF] Allibert | Estimating body-fixed frame velocity and attitude from inertial measurements for a quadrotor vehicle[END_REF], [START_REF] Leishman | Quadrotors and accelerometers: State estimation with an improved dynamic model[END_REF], [START_REF] Bangura | Non-linear velocity aided attitude estimation and velocity control for quadrotors[END_REF], [START_REF] Allibert | Velocity aided attitude estimation for aerial robotic vehicles using latent rotation scaling[END_REF].

If the drag coefficient for a given direction in the lateral plane is c > 0 and assuming no coupling drag terms between the lateral axes, then the drag force is given by [START_REF] Martin | The true role of accelerometer feedback in quadrotor control[END_REF], [START_REF] Mahony | Modelling and Control of Aerial Robots[END_REF] (see Figure 6)

D = - √ T   c 0 0 0 c 0 0 0 0   V.
If the accelerometer measurement is ā, and it measures the external forces i.e. the -T m e 3 + D m terms of (2a), then [START_REF] Leishman | Quadrotors and accelerometers: State estimation with an improved dynamic model[END_REF], [START_REF] Allibert | Estimating body-fixed frame velocity and attitude from inertial measurements for a quadrotor vehicle[END_REF], [START_REF] Bangura | Non-linear velocity aided attitude estimation and velocity control for quadrotors[END_REF]

ā = - 1 m -D T ,
where the model for the drag force in any direction is given by

D e 1 = -c √ T V ⊤ e 1 , D e 2 = -c √ T V ⊤ e 2 and D ⊤ e 3 = 0. Given that T = -mā ⊤ e 3 , implies āx = -c √ -ā z V x , āz < 0, (3a) āy = -c √ -ā z V y , āz < 0 (3b)
where c is now scaled by 1 √ m . These equations are used to obtain the translational linear velocity measurements in {B} [START_REF] Bangura | Non-linear velocity aided attitude estimation and velocity control for quadrotors[END_REF].

To obtain the vertical measured velocity B V z , an approach proposed in [START_REF] Allibert | Velocity aided attitude estimation for aerial robotic vehicles using latent rotation scaling[END_REF] which uses the aerodynamic power developed in [START_REF] Bangura | Aerodynamic power control for multirotor aerial vehicles[END_REF] could be used. Such an approach though very accurate requires measurement of electrical power and speed of all the rotors. In general, rotor speed measurement is absent on most generic open-source multirotor electronic speed controllers (ESC) though in recent years, there has been an increasing number of flight control boards that are equipped with voltage and current sensors for measuring electrical power. However, the majority of flight control boards are equipped with barometeric pressure sensors for measuring altitude. Using the altitude measurements z of the strapdown barometeric pressure sensors, we propose a complementary filter to obtain the measured vehicle velocity vz in e 3 in the inertial-fixed frame Filter 0.

ż = vz -k 1 (ẑ -z) , (4a)
vz = e ⊤ 3 Rā + g -β az -k 2 (ẑ -z) , (4b)
βaz = k 3 (ẑ -z) , ( 4c
)
where ẑ(0) = 0, k 1 , k 2 , k 3 are positive scalar gains, β az ∈ R is an estimate of the accelerometer bias in e 3 and ẑ is the estimated altitude.

Using the B Vx , B Vy obtained from ( 3) and ( 4), then using simple rotations based on the estimated attitude R, the full body-fixed frame velocity measurements B V can be obtained. Note that we assume for Filter 0 that the estimated attitude is correct R = R, however, small errors in attitude will not generate significant errors in vertical velocity. It should be noted that the concept of using B Vx , B Vy , A vz has been used by Hua et. al. [START_REF] Hua | Riccati nonlinear observer for velocity-aided attitude estimation of accelerated vehicles using coupled velocity measurements[END_REF] in their proposed filter along with proof of stability of the filter.

C. Inertial reference frame measurements of linear velocity

A vOur definition of an inertial-fixed frame sensor is a sensor that provides measurement such that it is independent of the relative airflow around the vehicle. For this work, the inertial-fixed frame linear velocity measurements are provided by either a uBlox global positioning system (GPS) or a Vicon motion capture system 1 . However, there are many other sensors and algorithms that can be used to obtain inertial reference frame measurements. Some of these algorithms, their required sensors and characteristics are surveyed in Table I. With measurements of linear velocity in both the bodyfixed and inertial frames available the challenge is to fuse these and provide estimates of attitude and velocity to be used by the controller proposed in Section III-A which is agnostic to inertial velocity measurements. It should be noted that the inertial measurement unit (IMU) on multirotors can provide measurements of both linear acceleration ā and angular velocity Ω. For this work we assume that the quadrotor is equipped with a magnetometer which measures the earth's magnetic field μ ∈ R 3 ∈ {A}. If R is the attitude of the vehicle, then the decoupled estimated magnetic field is μ = R⊤ μ.

1 https://www.vicon.com/

D. Quadrotor state estimation

The architecture of the proposed coupled non-linear complementary state estimation scheme is presented in Figure 2. In the filter formulation, please note

R 1 = B R = A R = R 2 = R and that B V = V 1 , A v = v 2 .
Before presenting the filter mathematically, we first define the anti-symmetric projection operator in square matrix space as

P a (H) = 1 2 (H -H ⊤ ). (5)
If we also let β ∈ R 3 to denote the bias in a measurements, subscripts 1 and 2 denote internal variables used by Filter 1 in {B} and Filter 2 in {A} respectively. If we assume that the trajectory of the quadrotor is smooth with Ω(t), Ω(t), v(t), v(t), V (t) and V (t) and bounded with no delays in the measurements, we propose the following coupled non-linear velocity aided attitude filters in both the body-fixed and inertial-fixed frames. Starting with the body-fixed frame {B} (Filter 1).

Filter 1 in {B}. Let k vc 1 , k r 1 , k u 1 , k βa
1 be positive scalar gains with measurements Ω = Ω, a = ā, V = V and µ = μ, then consider the following observer based on (2)

V1 = -(Ω -βΩ 1 ) × V1 + gX ⊤ 1 e 3 + a -βa 1 -∆ v 1 -k vc 1 ∆ vc 1 , (6a)
Ẋ1 = u1 R1 + u 1 Ṙ1 , (6b) Ṙ1 = R1 (Ω -βΩ 1 ) × -k r 1 ∆ r 1× R1 + R1 ∆ µ 1 -∆ rc 1 R1 , (6c) u1 = -k u 1 g Ṽ ⊤ 1 R⊤ 1 e 3 , ( 6d
) βa 1 = k βa 1 ∆ v 1 -k βac 1 ( βa 1 -βa 2 ), (6e) βΩ 1 = k βΩ 1 ∆ r 1 , (6f)
where Ṽ = V -V , X = u R is the scaled rotation with u a positive scalar. The initial conditions are

V1 (0) = V 1 (0), X 1 (0) = u 1 (0) R1 (0), u 1 (0) = 1, R1 (0) = I and I is the identity matrix. The innovation terms ∆ v 1 , ∆ r 1 , ∆ µ 1 are defined by ∆ v 1 = k v 1 ( V1 -V 1 ), (7a)
∆ r 1 = g u 1 R1 Ṽ1 × e 3 , (7b)
∆ µ 1 = k µ 1 (((µ × μ) ⊤ X ⊤ 1 e 3 )X ⊤ 1 e 3 ) × , (7c)
and the coupling innovation terms ∆ vc 1 and ∆ rc 1 are defined by

∆ vc 1 = V1 -R⊤ 1 (v 2 -ŵ), (8a)
∆ rc 1 = k rc 1 P a ( R1 R⊤ 2 ), (8b)
then for almost all initial conditions,

V1 (t) = V 1 (t) and R1 (t) = R 1 (t).
The stability analysis of this body-fixed frame filter (Filter 1) without any of the coupling and bias terms can be found in the authors' previous paper [START_REF] Allibert | Velocity aided attitude estimation for aerial robotic vehicles using latent rotation scaling[END_REF].

For the inertial frame filter (Filter 2), the following is proposed.

= v. Let k r 2 , k u 2 , k w 2 , k βa 2 , k v
2 be positive scalar gains, then consider the following observer

v2 = R2 (a -βa 2 ) + g e 3 -∆ v 2 -∆ vc 2 , (9a)
Ẋ2 = u2 R2 + u 2 Ṙ2 , (9b) Ṙ2 = R2 (Ω -βΩ 2 ) × -k r 2 ∆ r 2× R2 + R2 ∆ µ 2 -∆ rc 2 R2 , (9c) u2 = -k u 2 gṽ ⊤ 2 e 3 , (9d)
ẇ = -k w 2 ( ŵ -v2 ) + R1 V1 , (9e)
βa 2 = k βa 2 ∆ v 2 -k βac 2 ( βa 2 -βa 1 ), (9f) βΩ 2 = k βΩ 2 ∆ r 2 , ( 9g
)
with initial conditions v2 (0) = v 2 (0), u 2 (0) = 1, ŵ(0) = w(0), R2 (0) = I and ṽ(t) = v(t) -v(t). The innovations ∆ v 2
and ∆ r 2 are given by

∆ r 2 = g u 2 ṽ2 × e 3 , (10a)
∆ v 2 = k v 2 (v 2 -v), (10b) ∆ µ 2 = k µ 2 (((µ × μ) ⊤ X ⊤ 2 e 3 )X ⊤ 2 e 3 ) × . ( 10c
)
The coupling innovation terms ∆ vc 2 and ∆ rc 2 are given by

∆ vc 2 = k vc 2 (v 2 -ŵ) -R1 V1 , (11a)
∆ rc 2 = k c 2 P a R2 R⊤ 1 , (11b)
then for almost all initial conditions, v2 (t) → v 2 (t), ŵ → w and R2 (t) → R 2 (t).

The outputs of the coupled filter are the attitude estimates R, linear velocities B V , A v in the body-fixed and inertial frames of the vehicle respectively and wind velocity ŵ ∈ R 3 expressed in {A} at the 200Hz sampling frequency of the IMU. In addition, the filters also output the latent rotation scaling u for each filter. Given that the optimal value of u is unity, it is a measure of the performance of the filter in terms of minimising the errors between measurements and estimates. Though we do not provide stability proof as well as robustness of the filter, the validation is the results obtained in Section V. Also the positive scalar gains

k vc 1 , k r 1 , k u 1 , k βa 1 and k r 2 , k u 2 , k w 2 , k βa 2 , k v
2 are tuned based on classical pole placement technique on the linearised closed-loop system.

Remark 1. In certain situations such as indoors where the high ferrous content affects the magnetometer readings and thus the innovation ∆ µ 1 and ∆ µ 2 which use these readings can lead to errors in R e 3 . In such situations, it is necessary to set

k µ 1 = k µ 2 = 0.
It should be noted that due to the decoupling of our filter formulation through Equation 7c and 10c, these effects do not affect e 1 , e 2 and thus the overall stability of the filter.

E. Absence of inertial reference frame linear velocity measurements

In the absence of inertial reference velocity measurements v2 , the wind estimate ŵ is no longer observable. We propose to continue to implement (9e), however, we remove both the driving term v2 (which is not available) as well as the R1 V1 term to give ẇ = -k w 2 ŵ.

This ensures that there is no step change in the wind estimate ŵ and allows the estimate to decay to zero at the same rate at which the wind estimate is computed. If the wind estimate is explicitly used in perceptual control algorithms such as teleoperation then this approach ensures that the pilot has the opportunity to compensate and adjust accordingly. The proposed approach can also be justified in the sense that the normal scenario where GPS velocity information is lost corresponds to the case where the vehicle enters an enclosed space and the constant wind disturbance drops to zero in the new space.

III. SUPERVISORY CONTROL

In this section, we use the estimated velocities from Section II to do supervisory control of the vehicle. We propose a velocity controller in the body-fixed frame that is locally exponentially stable. From Filter 1 and Filter 2, we showed that estimates of the body-fixed frame linear velocities are always available irrespective of whether inertial reference frame velocity measurements are available or not. With this, the vehicle can be controlled in a supervised manner in {B} which is the frame in which humans perceive motion. We also show that if there is a need for an inertial reference frame control of the vehicle, the input to the proposed controller should be computed based on the estimated velocity of the wind.

A. Vehicle control system

The proposed (quadrotor) control system in the body-fixed frame {B} is shown in Figure 3. It has a hierarchical control structure with: high level velocity, middle level attitude and low level thrust control.

Consider the quadrotor dynamics in {B} in freestream with no forces acting on the vehicle given by (2) where V = B V and let the desired velocity set point in {B} be V d ∈ R 3 provided by the master haptic device. The resulting error Ṽ = V -V d has dynamics given by

V = -Ω × V + gR ⊤ e 3 - T m e 3 - T m cV + δ -Vd . ( 12
)
The following velocity controller which ensures local exponential stability of the resulting velocity error dynamics is proposed.

Controller 1. Let K i , K v ∈ R 3×3 be
mgR ⊤ d e 3 = Td e 3 + T d cV d + mΩ × V d -mK v Ṽ -m δ + m Vd , (13a) δ = K i Ṽ , δ(0) = δ0 , ( 13b
)
where δ ∈ R 3 is bounded and is an estimate of the model error δ and δ0 is some initial condition. The error dynamics of (12) are globally asymptotically and locally exponentially stable around the equilibrium Ṽ = 0 under the assumption that R ≈ R d and T ≈ T d .

The exponential stability proof of this theorem can be found in the authors' previous work [START_REF] Bangura | Non-linear velocity aided attitude estimation and velocity control for quadrotors[END_REF]. The desired attitude R d and thrust T d are passed as set points to a computationally efficient quaternion tracking controller [START_REF] Bangura | An open-source implementation of a unit quaternion based attitude and trajectory tracking for quadrotors[END_REF] acting as a middle level controller. As long as the gains of this controller are well tuned (using pole placement of the closed loop dynamics), the tracking errors are small even for aggressive motion demands due to the low rotational inertia and high actuation levels on the quadrotor vehicle. In practice, the acceleration Vd is not available and we set Vd = 0 without observing any notable performance changes. The limits of the controller are defined by the 2ms PWM signal sent to each speed controller and the rise time of the rotors.

In addition, if an outside observer in the inertial frame wants to control the vehicle, then the control should be done in the inertial frame. However, the vehicle should still be controlled in {B}, with the user input A v d written in {B} i.e. B V d using the following relationship

B V d = R⊤ A v d -ŵ .
This is then the input to the velocity controller Controller 1. This is very important since the assumptions in the models in Equations 1 and 2 can be broken (e.g. when the freestream assumption is broken) resulting in a wrong B V . However choosing to control using an outside observer reference, the wrong estimate of w and B V results in the correct A v.

Remark 2. Given that the proposed architecture is for any generic multirotor electronic speed controller, in order to minimise the controller dependence of III-A on δ e 3 error term, the following PWM to battery voltage scaling is used

T c = V O -Vb V O T d ,
where V O is some nominal battery voltage that is a function of the number of cells of the battery, Vb is the instantaneous measured bus voltage, T d is the computed desired thrust and T c is the actual commanded thrust as PWM signal sent to the speed controller.

IV. HARDWARE SETUP In order to do the supervisory control, we choose to use force feedback teleoperation as it addresses situations where high bandwidth feedback of velocity is required for dynamically controlling vehicles that dynamically react with the environment. Furthermore, teleoperation generally does not take into account feedforward information of a trajectory and is applicable to situations where position measurements may not be available. This section describes the hardware setup for the experimental results presented in Section V. The setup comprises of the following components: master joystick device, ground station, communication link and multirotor robot. In addition, a radio controller is used for switching between manual and autonomous control as a safety precaution. The entire hardware setup is shown in Figure 4.

A. Master haptic joystick for supervisory control

An admittance configured bilateral force/haptic feedback teleoperation scheme is employed in this paper as the master device. The master device is used to measure force inputs from the pilot and exert motion feedback. Our approach was tested with a custom built admittance joystick shown in Figure 5(a in Figure 5(b). This commercial off-the-shelf joystick ensures that we maintain generality and the results of this paper can be replicated. The NOVINT Falcon employs a delta robot mechanical structure that is capable of measuring the position of the joystick's end effector and exerting a scaled force. A high gain PID (Proportional Integral Derivative) position controller and force estimator are implemented to approximate the admittance configuration [START_REF] Alaimo | Admittance-based bilateral teleoperation with time delay for an unmanned aerial vehicle involved in an obstacle avoidance task[END_REF], [START_REF] Hou | Representation of vehicle dynamics in haptic teleoperation of aerial robots[END_REF]. Given the limitations of the joystick hardware, the translational degrees of freedom of the joystick are mapped to the linear velocities of the slave vehicle and the side buttons are used to change the heading of the robot.

B. Ground station

The master haptic device controller, Vicon data routing module, control mode module and logging modules are implemented on the ground station using ROS (Robot Operating System) 3 . A pair of 915MHz 3DR radio transceivers are used for communication between the ground station and the robot at 40Hz using MAVLINK protocol.

C. Slave robot

A custom built quadrotor shown in Figure 1 is used as the slave robot. It uses the Pixhawk flight control board which has an MPU6000 IMU that provides measurement of linear acceleration ā and vehicle angular velocity Ω, a HMC5883 magnetometer for measuring the earth's magnetic field μ [START_REF] Meier | Px4: A node-based multithreaded open source robotics framework for deeply embedded platforms[END_REF] and a MEAS MS5611 barometer. 3 http://www.ros.org/

V. EXPERIMENTAL RESULTS

In this section, we present the experimental results of the proposed scheme. The section starts by presenting the calibration results for the body-fixed frame velocity measurements. These results (both attitude and linear velocities) are compared to Vicon motion capture ground truth measurements. The state estimation results are then presented for low and high velocities as well as the teleoperation results for flying the slave multirotor vehicle with and without inertial reference frame measurements. We also present results of wind velocity estimation.

It should be noted that for the drag force to IMU calibration, in order to minimise notational confusion, the accelerometer measurements used here are ā = āβ a 1 , where β a ∈ R 3 was defined in Section II-D as the bias estimated from the attitude observer in Filter 1. Given also the slow dynamics of the vehicle compared to the filter, we make the following necessary assumption βa 1 = β a 1 . āβ a 1 is as close as one can get to the true acceleration of the vehicle. This is what distinguishes our work from previous work and gives us better results than the state-of-the art [START_REF] Burri | Robust state estimation for micro aerial vehicles based on system dynamics[END_REF], [START_REF] Sikkel | A novel online model-based wind estimation approach for quadrotor micro air vehicles using low cost mems imus[END_REF], [START_REF] Davis | Direct sensing of thrust and velocity for a quadrotor rotor array[END_REF], [START_REF] Tomić | The flying anemometer: Unified estimation of wind velocity from aerodynamic power and wrenches[END_REF], [START_REF] Arain | Real-time wind speed estimation and compensation for improved flight[END_REF], [START_REF] Leishman | Quadrotors and accelerometers: State estimation with an improved dynamic model[END_REF].

A. Drag calibration for B V translational velocities

The major calibration scheme required is for the drag coefficient c. Generally there is an offset in the accelerometer measurement due to mass imbalance and from each power cycle. Hence the drag force model (3), can be written to incorporate this offset say for the e 1 direction as

āx -a x of f = -c -ā z + a z of f V x ,
where a x of f and a z of f are the offsets in the accelerometer measurements. For ease of applying linear regression, this model is simplified to

āx √ -ā z = -(cV x + ǭ) , āz < 0,
where ǭ is the new offset. This gives rise to two types of calibration for the drag coefficient: offline and online. In the offline process, both c and ǭ are determined while the online case involves recalibrating to determine the mass imbalance and the new offset ǭ that is as a result of β a which occurs on each power cycle. This offset ǭ is determined when inertial reference frame linear velocity measurements are available and at the start of every flight when inertial reference linear velocity measurement is available. This is done by commanding the vehicle to hover i.e. V d = 0. Furthermore using GPS in an outdoor environment, we tested the first order drag model along the e 1 axis to show that the linear model for accelerometer measurement to velocity holds even at speeds of up to 14m/s. The result is presented in Figure 6. The figure also has a comparison between the linear drag model and the quadratic model used to capture parasitic drag forces. The resulting regression coefficients for both the linear and quadratic models are

āx = - √ -ā z 0.0794 Vx + 0.0224 , āx = - √ -ā z 0.00209 Vx | Vx | + 0.02586 Vx + 0.0162 , -5 0 5 10 V x [m/s] -1 0 1 2 ax √ -az

Drag coefficient determination

Linear R 2 = 0.4272 Square R 2 = 0.47436

Vx Vx 2 Vx
Fig. 6. Drag coefficient c determination using linear and quadratic drag models. The blue represent the raw GPS measurements, red line is the linear model and black is the quadratic model.

with R 2 goodness of fit 0.4272 and 0.4736 respectively. Hence even at speeds of up to 14m/s, the linear drag force model holds as the contribution of the parasitic drag term 0.00209 V 2

x in the quadratic model is small compared to the other drag terms described in [START_REF] Bangura | Nonlinear dynamic modeling for high performance control of a quadrotor[END_REF]. These drag terms were lumped into a linear model which we have defined to be a function of c. Hence even at higher velocities, the dominant drag forces are flapping and induced drag. It is important to note that the low R 2 goodness of fit values for Figure 6 are as a result of noisy GPS inertial reference frame velocity measurements rotated into the body-fixed frame. This linear modelled is further verified by the high speed results obtained using GPS shown in Figure 8. Using the obtained regression coefficients for the drag constant and offsets the measurement of planar velocities of the vehicle in the body-fixed frame are obtained. With these measurements, the estimated attitude of the vehicle, and using barometer measurements along with Filter 0, the full bodyfixed frame measured velocities is obtained. Figure 7 shows the resulting computed velocities (red) and ground truth Vicon motion capture measurements (blue). From these results, one can conclude that the computed velocities match those of the ground truth and thus can be used as measurements for bodyfixed frame linear velocities B V .

B. State estimation results

To show the effectiveness and accuracy of the estimated attitude and velocities, we compare our estimates to Vicon motion capture ground truth measurements. A comparison of the attitude estimates to the ground truth attitude is shown in Figure 9(a). From this, it is clear that there is negligible attitude error between the ground truth and the estimated attitude. It should be noted that the latent scaling u though its optimal value is 1, the time variations are as a result of errors between the measurements and estimated states. Given that the velocity measurements in the body-fixed frame came in at 200Hz and with significantly very high noise implies that if one is to use the estimates in a controller, obtaining smooth velocity estimates are more important than minimising the total error in the velocity estimates. The resulting estimates for velocities State estimation results for high speed using GPS is shown in Figure 8. A summary of the mean error and standard deviation of the errors is presented in Table IV. To show superiority of our filter, Table III compares the errors of the outputs from the proposed filter to the current state-of-the art.

To prove that our scheme can estimate wind velocity when inertial measurements are available, a set of four fans with streamers were placed (x ≈ 3.0) at the end opposite the vehicle's starting location and were turned on while the vehicle was flown with closed loop control in e 1 and e 2 directions while the e 3 was manually controlled as shown in Figure 10. The flow field generated by the fans is highly turbulent and 11 and mean and covariance summarised in Table IV. Though the flow was highly turbulent, our best attempt at obtaining the average velocity measurements using a Mastech ms6252b digital anemometer 4 in the longitudinal direction is also shown along with our estimated results. From these results, we can conclude that our scheme is able to estimate the wind velocity in e 1 , e 2 directions. It should be noted that the high turbulent nature of the flow requires a gust model which is beyond the scope of this work. Furthermore, the computational efficiency expressed as a percentage of CPU usage of the proposed filter compared to other filters implemented on the PX4 stack firmware [START_REF] Meier | Px4: A node-based multithreaded open source robotics framework for deeply embedded platforms[END_REF] is shown in Table IV. It is clear that the SO(3) filter is more computationally efficient than our proposed filter as it is implemented in quaternions. However, this filter only outputs attitude. Running this along with the local position filter to provide estimates of velocity and attitude, results in a low computationally efficient filter combination than our proposed complementary filter. Furthermore, the extended Kalman filter (EKF) based attitude and position filter is twice more computationally expensive than the proposed filter.

C. Teleoperation with and without Inertial Reference measurements

To demonstrate the effectiveness of the scheme presented, we teleoperated the vehicle in a Vicon motion capture arena and turned off the Vicon motion measurements (from time (b) Ground truth and estimated velocities expressed in the body-fixed frame. Blue is true inertial reference measurement rotated into the body-fixed frame R ⊤A v. Red is estimated inertial rotated into the body-fixed frame using the filter rotation estimate R⊤A v. Green is estimated body-fixed frame velocity B V . In this experiment the wind velocity is w = 0.

Fig. 9. Comparison of attitude and velocity estimates from proposed complementary filter structure for low speeds using Vicon as the reference system. t = 146 and t = 200) to simulate loss of inertial reference measurements. The results of the scheme working with and without Vicon inertial reference measurements are shown in Figure 12 and 13. From Figure 13, the wind estimates, it is clear that the vehicle was teleoperated with and without (between t = 146 and t = 200) any inertial reference velocities. From these results, it is clear that the scheme works whether there is inertial reference velocity measurements or not and the entire observer-controller does not go unstable when inertial measurements become available in mid air. A video 5 of one of the flight experiments is available online. In this paper we have formally presented our proposed coupled non-linear complementary velocity aided attitude filter that ensures the availability of estimated body-fixed frame linear velocities under all circumstances. A drag force to accelerometer model along with strapdown barometer measurements are used to provide velocity measurements in the body-fixed frame while inertial reference frame sensors such as Vicon, GPS, vision and laser based Inertial Navigation Systems (INS) provide linear velocity measurements in the inertial reference frame. Results were presented to show comparisons to ground truth Vicon and GPS measurements. Results were also presented to demonstrate the capability of the scheme to estimate the wind velocity when inertial measurements are available. To test the proposed scheme, an admittance configured bilateral teleoperation system was used in supervisory control of the vehicle. This allowed a human pilot to provide velocity set points derived from force inputs. These desired velocities were then tracked by an onboard locally exponentially stabilising velocity controller. A performance comparison of the proposed filter to other filters and most recent techniques shows promising results in terms of computational efficiency and accuracy in estimating low and high linear velocities. As a future work, we will provide a full proof of the coupled filter and derive the quaternion version of the proposed filter which will result in a further increase in the computational efficiency of the scheme.

Fig. 1 .

 1 Fig. 1. Quadrotor platform: The platform used for the indoor experiments showing the body-fixed {B} and inertial {A} frames.

Fig. 2 .

 2 Fig. 2. The proposed coupled non-linear complementary velocity aided attitude filter. Input to the filter are inertial reference frame velocity, IMU, magnetometer and barometer. The coupled filters are Filter 1 in {B} and Filter 2 in {A}. The diagram also shows Filter 0 which uses barometer measurements to provide vertical body-fixed frame velocity measurements.

Filter 2 in

 2 {A}. Consider the model (1) in the inertial-fixed frame {A} with wind speed w and measurements a = ā, Ω = Ω and v

Fig. 3 .

 3 Fig. 3. Proposed hierarchical control structure on the vehicle.

Fig. 4 .

 4 Fig. 4. Hardware layout showing the various components: FALCON, ground station and slave robot.

Fig. 5 .

 5 Fig. 5. Master haptic devices used in this paper. A custom made admittance and a COTS NOVINT Falcon joystick.

Fig. 8 .

 8 Fig.8. Filter outputs compared to ground truth for high-speed outdoor flight in almost zero wind velocity. Ground truth and estimated velocities expressed in the body-fixed frame. Blue is true GPS inertial reference velocity rotated into the body-fixed frame R ⊤A v. Red is estimated inertial rotated into the body-fixed frame using the filter rotation estimate R⊤A v. Green is estimated body-fixed frame velocity B V .

Fig. 10 .

 10 Fig. 10. Quadrotor flying in synthetic wind field for wind velocity estimation. The red streamers show the direction and strength of the airflow.

Fig. 11 .

 11 Fig. 11. Results for estimated wind velocities (red) compared to measured average wind velocities (blue) obtained from the Mastech ms6252b digital anemometer. The experiments were performed for a total of 10 times with the average shown

  Velocity controller results, further demonstrating the accuracy of the estimated body-fixed frame velocities to Vicon measurements. In the figures, the blue lines indicate the desired velocity set point, red is the estimated inertial reference velocity in the body-fixed frame and black is the estimated body-fixed frame velocity using the drag model.

  Force input and position response of haptic joystick. Red is scaled user input force, blue is joystick output force.

Fig. 12 .

 12 Fig. 12. Teleoperation with and without Vicon inertial reference frame measurements. The Vicon measurements were turned off i.e. k w 2 = ŵ = 0 between t = 146 and t = 200.

Fig. 13 .

 13 Fig. 13. Estimated wind velocities during the flight with no synthetic additional wind. Between t = 146 and t = 200, the Vicon measurements were turned off.

TABLE I INERTIAL

 I REFERENCE FRAME ESTIMATION ALGORITHMS FOR GPS DENIED ENVIRONMENTS ALONG WITH THEIR LIMITATIONS. THOUGH THESE SENSORS USUALLY MEASURE IN THE CAMERA FRAME, THEY ACTUALLY DO SO WITH ZERO WIND, HENCE, CAN BE CONSIDERED INERTIAL FRAME SENSORS.

TABLE II COMPARISONS

 II OF FILTERS AVAILABLE ON THE PX4 STACK.

	Mode	CPU computa-	Features
		tional cost (µ ±
		σ)
	SO(3)	3.1062±0.043%. Quaternion version of the non-
			linear complementary SO(3) filter
			in [49].
	Attitude EKF	14.916±0.205%. Extended Kalman filter for atti-
			tude only estimation.
	Position Estima-	10.846±0.21%.	Extended Kalman filter for posi-
	tor		tion and velocity estimation.
	Attitude and Po-	22.988±0.374%. Extended Kalman filter for esti-
	sition EKF		mating position and attitude.
	Proposed filter	10.711	Non-linear complementary veloc-
		±0.075%.	ity aided attitude filter.

  Estimated attitude (blue) compared with Vicon ground truth (red). One should note that the errors in the latent scaling u are as a result of high noise of the velocity measurements in {B} which resulted in using low k v 1 gain.

						Estimated and true	φ
	[deg] φ	20 0 -20									φ φ
			65	70	75	80	85	90	95	100	105
							Time [s]
						Estimated and true	θ
	[deg] θ	-10 10 0									θ θ
			65	70	75	80	85	90	95	100	105
							Time [s]
						Estimated and true	ψ
	[deg] ψ	105 115 120 110									ψ ψ
			65	70	75	80	85	90	95	100	105
							Time [s]
						Latent scaling u
	u	0.8 1.2 1									u d û
			65	70	75	80	85	90	95	100	105
							Time [s]
	-2 2 (a) 65 V 0 x [m/s]	70	75	80	85	90	95	100	105 R ⊤A vx R⊤A vx B Vx
							Time [s]
		2
	y [m/s]	0									R ⊤A vy R⊤A vy B Vy
	V
		-2	65	70	75	80	85	90	95	100	105
							Time [s]
		2
	z [m/s]	0									R ⊤A vz R⊤A vz B Vz
	V
		-2	65	70	75	80	85	90	95	100	105
							Time [s]

4

http://www.p-mastech.com/product/detail/426

TABLE IV VELOCITY

 IV ERROR STATISTICS IN M/S OF OUR PROPOSED FILTER

	Variable	µerror	σerror
	wx	0.1997	0.0724
	wy	0.0320	0.0107
	wz	0.0434	0.0021
	Low ṽx 1.6m/s	0.1349	0.0210
	Low ṽy 1.6m/s	0.1291	0.0281
	Low ṽz 1m/s	0.1552	0.0428
	High ṽx 14.8m/s	1.3713	1.1818
	High ṽy 8m/s	0.7676	0.9676
	High ṽz 5m/s	1.6571	1.9896

https://www.youtube.com/watch?v= ntGvVv7eyw

goo.gl/KNsuXE


• Objective : characterize the impact of the "inexact" mode-splitting on the stability of numerical models

Our approach :

1. Stability analysis of the mode splitting technique based on an eigenvector decomposition using the exact barotropic mode 2. Quantify the amount of dissipation required to stabilize the approximative splitting 2D linearized primitive equations (u 0 = w 0 = 0, f cor = 0)

                   ∂tu + 1 ρ0 ∂x p = 0 ∂z p = -ρg ∂xu + ∂zw = 0 ∂t ρ + w dρ dz = 0
with boundary conditions

w(z = 0) = ∂tη η H w(z = -H) = 0 (flat bottom) p(z = 0) = ρ0gη F. Lemari é -Recent Developments around CROCO 2D linearized primitive equations (u 0 = w 0 = 0, f cor = 0)                    ∂tu + 1 ρ0 ∂x p = 0 ∂z p = -ρg ∂xu + ∂zw = 0 ∂t ρ + w dρ dz = 0
with boundary conditions

w(z = 0) = ∂tη η H w(z = -H) = 0 (flat bottom) p(z = 0) = ρ0gη Sturm-Liouville eigenvalue problem with Λ = - d dz N -2 d dz            ΛMq(z) = c -2 q Mq(z) dMq dz (-H) = 0 dMq dz (0) = - N 2 (0) g Mq(0) Mq, Mm = δq,m w.r.t. the inner product f, g = 1 H 0 -H f gdz F. Lemari é -Recent Developments around CROCO 2D linearized primitive equations (u 0 = w 0 = 0, f cor = 0)                    ∂tu + 1 ρ0 ∂x p = 0 ∂z p = -ρg ∂xu + ∂zw = 0 ∂t ρ + w dρ dz = 0
with boundary conditions

w(z = 0) = ∂tη η H w(z = -H) = 0 (flat bottom) p(z = 0) = ρ0gη
Expand the variables u, p in the eigenfunctions M q (z) of

Λ              u(x, z, t) = K-1 q=0 uq(x, t)Mq(z) p(x, z, t) = gρ 0 K-1 q=0 hq(x, t)Mq(z)
to obtain a set of uncoupled systems

     ∂tuq + g∂xhq = 0 ∂thq + c 2 q g ∂xuq = 0
In particular

• q = 0 → barotropic mode

• q > 1 → baroclinic modes Starting from the PEs, horizontal velocities are decomposed as

u = 1 H η -H u dz + u 3d = u + u 3d , η -H u 3d dz = 0
Standard barotropic mode (associated to the 2D linearized PEs)

Integrating vertically the continuity and momentum equations we end up with

   ∂tη + ∂x(Hu) = 0 ∂tu + g∂xη = - 1 ρ0 ∂x 1 H 0 -H p h dz
Remark : the red term is a "slow" term absent from the normal mode analysis.

→ In practice it is kept frozen during the barotropic integration

F. Lemari é -Recent Developments around CROCO

Interpretation in terms of normal modes

The underlying assumption is that the barotropic mode is depth-independent

M 0 (z) = 1, u 0 (x, t) = u, 0 -H M q (z)dz = 0 (q ≥ 1)
The orthogonality condition between modes is lost → Exact normal mode decomposition → K independent systems → Depth-independent barotropic mode assumption → K -1 independent systems (for q ≥ 1) + a system which includes contributions from all modes F. Lemari é -Recent Developments around CROCO Consequences of the depth-independent assumption 1. Some (fast) barotropic contributions are treated as "slow" terms

   ∂tη + ∂x(Hu) = 0 ∂tu + g∂xη = - 1 ρ0 ∂x 1 H 0 -H p h dz 1 H 0 -H p h dz = 1 H 0 -H (p-ρ0gη) dz = ρ0g K-1 q=0 1 H 0 -H Mq(z) -Mq(0) hq
⇒ the contribution of the barotropic mode to the red term is proportional to

For N 2 = cste 1 H 0 -H M0(z) -M0(0) ≈ ε 3 , ε = N 2 H g 2.
The consistency between the barotropic and baroclinic modes must be enforced since they are no longer independent

Time stepping the coupled baroclinic-barotropic system

→ Dissipation is necessary to stabilize the time integration procedure because some fast contributions are treated with the slow time-step.

Objective : determine the minimum amount of dissipation required to stabilize the coupled baroclinic-barotropic system integration 1. Averaging filters (e.g. Nadiga et al., 1997;Shchepetkin & McWilliams, 2005)

u n+1,f = N filter m=0 amu m (N filter ≥ N split )
under the constraints

N filter m=0 am = 1, N filter m=0 am m N split = 1
Several options:

-Flat filter

-Cosine filter n n + 1 n + 2 n + 3/2 n + 1/2 am Flat filter

Cosine filter

With averaging filters -the amount of dissipation is somewhat ad-hoc (no quantification of its impact) -the barotropic integration must go beyond n + 1 -the barotropic components are not continuous in time F. Lemari é -Recent Developments around CROCO 2. Dissipative time-stepping (e.g. Hallberg, 1997) In the case of the CROCO ocean model : → a generalized forward-backward scheme (AB3-like, AM4-like) is used for the barotropic mode (on top of an averaging filter)

η n+1 = η n -H∆t∂x (3/2 + β)u n -(1/2 + 2β)u n-1 + βu n-2 u n+1 = u n -g∆t∂x δη n+1 + (1 -δ -γ -)η n + γη n-1 + η n-2
-Revised choice : no averaging filter + additional constraint for the choice of (δ, γ, β, ) to guarantee stability of the split-explicit formulation

H(x) ⌘(x, t = 0) ⇢(x, t = 0) N = 1.4 ⇥ 10 3 s 1 x
In the nonlinear case there are other sources of instability :

• The rhs of the barotropic equations is integrated from -H to η

• Aliasing errors when the advective terms are added

F. Lemari é -Recent Developments around CROCO
Barotropic mechanical energy

E(t) = L 0 ρ0 2 (H + η)u 2 dx + L 0 ρ0 2 gη 2 dx
Comparison of standard ROMS vs modified CROCO scheme (for a fixed baroclinic time-step) E(t) / E(0)

p = pa + pH + q, pH = ρ0gη + g 0 z (ρ bq -ρ0)dz
Boundary conditions :

∂tη = w(0), w(-H) = 0, p(z = 0) = ρ0gη → q(z = 0) = 0
q : non-hydrostatic pressure which cancels the divergent part of velocity field

Implicit algorithm

1. Advance η with q = 0 or q

= q n η n+1 -η n ∆t -gH∆t∂ 2 x η n+1 = -H∂xu n + ∆tH ρ0 ∂ 2 x q
2. Compute provisional velocity field

u n+1 = u n -∆tg∂xη n+1 , w n+1 = w n 3. Solve ∆q = ρ0 ∆t ∂x u n+1 + ∂z w n+1
4. Correct velocity field to remove divergent part

u n+1 = u n+1 -∆t∂xq, w n+1 = w n+1 -∆t∂zq
But in the end :

η n+1 -η n ∆t = -H∂xu n+1
→ constancy preservation for tracers is lost + 2D/3D inconsistencies Implicit algorithm (Casulli; 1999)

1. Compute provisional η with q = 0 or q = q n η n+1 -η n ∆t -gH∆t∂ 2 x η n+1 = -H∂xu n + ∆tH ρ0 ∂ 2 x q . . .

5.

Correct free-surface to satisfy B.C.

η n+1 -η n ∆t = -H∂xu n+1
Still 2D/3D inconsistencies (may require a decrease of ∆t3D)

Homogeneous linearized equations

∂xu + ∂zw = 0 ∂tu = -g∂xη -∂xq/ρ0 ∂tw = -∂zq/ρ0 ∂tη = w(0) = -H∂xu Split-explicit algorithm 0 ≤ m ≤ N split -1 1.
Advance η and u with q = 0 or q = q n u m+1 = u m -g(δt)∂xη mδt ρ 0 ∂xq η m+1 = η m -δtH∂xu m+1 , 2. Compute provisional fields u n+1 and w n+1 5. Correct velocity field to remove divergent part

u n+1 = u n+1 -∆t∂xq, w n+1 = w n+1 -∆t∂zq However : u n+1 = u n+1
Homogeneous linearized equations

∂xu + ∂zw = 0 ∂tu = -g∂xη -∂xq/ρ0 ∂tw = -∂zq/ρ0 ∂tη = w(0) = -H∂xu Split-explicit algorithm 0 ≤ m ≤ N split -1 1.
Advance η and u with q = 0 or q = q n u m+1 = u m -g(δt)∂xη mδt ρ 0 ∂xq η m+1 = η m -δtH∂xu m+1 , 2. Compute provisional fields u n+1 and w n+1 5. Correct velocity field to remove divergent part

u n+1 = u n+1 -∆t∂xq, w n+1 = w n+1 -∆t∂zq However : u n+1 = u n+1
Solution 1 : change boundary condition on q to ∂zq| z=0 = 0

⇒ u n+1 = u n+1 = u n+1
Homogeneous linearized equations

∂xu + ∂zw = 0 ∂tu = -g∂xη -∂xq/ρ0 ∂tw = -∂zq/ρ0 ∂tη = w(0) = -H∂xu Split-explicit algorithm 0 ≤ m ≤ N split -1 1.
Advance η and u with q = 0 or q = q n u m+1 = u m -g(δt)∂xη mδt ρ 0 ∂xq η m+1 = η m -δtH∂xu m+1 , 2. Compute provisional fields u n+1 and w n+1 Semi-implicit forward-backward

u m+1 = u m -δt g∂xη m + c 2 s ∂xδρ m w m+1 = w m -δtc 2 s ∂z δρ m+θ δρ m+1 = δρ m -ρ0δt ∂xu m+1 + ∂zw m+θ η m+1 = η m + δt( w| z=0 ) m+θ
In practice :

• cs is chosen such that: ε = √ gH cs 1

• The acoustic waves are integrated, in a split-explicit free surface approach, at the same level (i.e. with the same time step) than the barotropic mode.

Advantages of the pseudo-compressible approach (Auclair et al., 2017)

• Surface waves "feel" the NH effects

• Main limiting CFL condition : horizontal acoustic waves propagation → may require to decrease the barotropic time-step

• This approach allows internal waves with high order divergence -pressure gradient computation

• Scales well with the resolution

• Compressibility effects physically significant ?

F. Lemari é -Recent Developments around CROCO

From nesting to truly multiresolution

What is needed for multiresolution ?

• Neighboring/overlapping grids connection

• Bathymetry smoothing

•F→

  numerical kernel Online nesting capability via Agrif library Non-hydrostatic Non-Boussinesq from S-NH Sediment module from Mars3d OAW coupling interface (shared w. Nemo) ALE-type vertical coordinate (ongoing + shared w. Nemo) Supported by Ifremer, Shom, Cnrs, Ird, Inria• Croco as a way to sustain the COMODO group Complementary analysis of the mode-splitting and control of numerical dissipation 2. A non-hydrostatic non-Boussinesq algorithm, the how and why 3. Multiresolution simulations using block structured mesh refinement 4. Future perspectives and concluding remarks Spurious "numerical mixing" not only associated with the space/time discretization of the advection operator

F

⇒Figure :

 : Figure : Time evolution of the maximum absolute value of the free surface elevation. Comparaison of usual filters against a reference solution without splitting.

  hydrostatic (Boussinesq) equations ∂tu = -∂xp/ρ0 ∂tw = -(∂zp -ρg) /ρ0 ∂xu + ∂zw = 0 ρ = ρ bq (θ, S, -ρ0gz) Pressure decomposition :

F

F

  ∂xu + ∂zw = 0 ∂tu = -g∂xη -∂xq/ρ0 ∂tw = -∂zq/ρ0 ∂tη = w(0) = -H∂xu

3.

  Correct u n+1 to enforce u n+1 = u n+1 4. Solve ∆q = ρ0 ∆t ∂x u n+1 + ∂z w n+1

3.

  Correct u n+1 to enforce u n+1 = u n+1 4. Solve ∆q = ρ0 ∆t ∂x u n+1 + ∂z w n+1

3.

  Correct u n+1 to enforce u n+1 = u n+1 4. Solve ∆q = ρ0 ∆t ∂x u n+1 + ∂z w n+1 5. Correct velocity field to remove divergent part u n+1 = u n+1 -∆t∂xq, w n+1 = w n+1 -∆t∂zq However : u n+1 = u n+1 Solution 2 : reset u n+1 = u n+1 and η n+1 -η n ∆t = -H∂xu n+1 ⇒ loose the gain in accuracy of the split-explicit algorithm Pressure correction / Pressure projection method • In generalized vertical coordinates, the Poisson pressure equation is hard to solve (25 diagonals, non symmetric) • Mode splitting: choices -relax the compatibility of 2D/3D fluxes (large source of instabilities) -relax the true boundary condition q(z = 0) = 0 -relax the accuracy of the 2D integration (for split-explicit algorithms) • Likely to require to decrease the baroclinic time-step compared to hydrostatic simulations 2D linearized NH Non-Boussinesq equations Taylor expansion of density field (with ∂ρ ∂p = c -2 s ) ρ = ρ(θ, S, p) = ρ bq (θ, S, p ref ) + ∂ρ ∂p δp δρ +O(δp 2 ) ∂tu = -∂xp/ρ0 ∂tw = -(∂zp -ρg) /ρ0 ∂tδρ = -ρ0(∂xu + ∂zw) Pressure decomposition : p = pa + pH + c 2 s δρ, pH = ρ0gη + g 0 z (ρ bq -ρ0)dz Boundary conditions : ∂tη = w(0), w(-H) = 0, p(z = 0) = ρ0gη → δρ(z = 0) = 0 Homogeneous linearized equations ∂tu = -g∂xη -c 2 s ∂xδρ ∂tw = -c 2 s ∂zδρ ∂tδρ = -ρ0(∂xu + ∂zw) ∂tη = w| z=0 w| z=-H = 0 δρ| z=0 = 0

-•FF1

  e.g. distance to the coast • Load balancing on parallel computers • Inputs / Outputs and visualization -Minimum number of grid points per grids : 256 (16x16) • Percentage of land points at the finer level : 20% • Cost of intergrid operations : 15% Sea surface temperature (multiresolution from ∆x = 800 m to ∆x = 100 m) percentage of execution time spent in cache memory accesses 2 floating point arithmetic instructions per Memory Read or Write → e.g. results used as a basis for the NEMO development strategy

F. Lemari é -Recent Developments around CROCO

Summary & perspectives on mode splitting analysis

The barotropic mode is traditionally considered as depth independent which is an assumption (e.g. work of Higdon, Bennett, de Szoeke) CROCO approach :

→ Provide a general framework for the stability analysis of the mode splitting approach → The framework allow the design of efficient 2D/3D time stepping algorithms → Approach implemented in standard Croco oceanic model Perspectives :

→ Transfer toward operational centers (in progress with Mercator-Ocean) Demange J., L. Debreu, F. Lemari é, P. Marchesiello, E. Blayo: Stability analysis of split-explicit free surface ocean models: implication of the depth-independent barotropic mode approximation, submitted to J. Comp. Phys.

⇒ Robust enough to be used for arbitrary realistic simulations

Feedback from model intercomparison Experience feedback from COMODO project

Positive aspects

• Inter-disciplinary collaborations around specific scientific questions

• Useful for code debugging and checking of numerical implementation

• Helpful to prioritize future code developments (e.g. Ilicak et al., 2011)

• Dynamic effect on the community with experience-sharing Difficulties/limitations • Difficult to define reference solutions

• "Good" solutions obtained for the wrong reasons

• Focused on dynamical core, no parameterization

• Generally simple geometries


Zuqi Tang

Patrick Dular

Yvonnick Le Ménach

Emmanuel Creusé

Francis Piriou

Comparison of Residual and Hierarchical Finite Element Error Estimators in Eddy Current Problems

Abstract⎯ The finite element computation of eddy current problems gives numerical error. This error cannot be calculated, but can only be estimated. Among all error estimators already developed, it is proposed to compare two proven estimators called residual and hierarchical error estimators.

I. INTRODUCTION

To evaluate the quality of numerical computation of low frequency electromagnetic finite element (FE) problems, a lot of error estimators have been developed. Some of them are based on energy minimization. Others, called equilibrated estimator, are based on the verification of the constitutive relations. More recently the residual a posteriori error estimator (REE) [START_REF] Babuska | A posteriori error estimates for the finite element method[END_REF][START_REF] Beck | Residual based a posteriori error estimators for eddy current computation[END_REF][START_REF] Chen | An adaptive finite element method forthe eddy current model with circuit/field couplings[END_REF][START_REF] Creusé | Residual-based a posteriori estimators for the A/Phi magnetodynamic harmonic formulation of the Maxwell system[END_REF] and the hierarchical error estimator (HEE) have been proposed [START_REF] Dular | A posteriori error estimation of finite element solutions via the direct use of higher order hierarchal test functions[END_REF]. In the case of the REE, all the terms weakly verified are evaluated to obtain an estimation of the error [START_REF] Creusé | Residual-based a posteriori estimators for the A/Phi magnetodynamic harmonic formulation of the Maxwell system[END_REF]. On the other hand, with the HEE, the error estimation is given by using hierarchical higher order (HO) test functions (TFs) [START_REF] Webb | Hierarchal scalar and vector tetrahedra[END_REF] in the FE formulation, either to evaluate some residues or to obtain HO solutions [START_REF] Dular | A posteriori error estimation of finite element solutions via the direct use of higher order hierarchal test functions[END_REF].

In this paper, we propose to apply the REE to FEs up to second order and to compare the obtained results with the HEE on the same mesh.

II. ERROR ESTIMATORS

A. Weak finite element magnetodynamic formulation

The weak magnetic vector potential a-formulation of the magnetodynamic problem is obtained from the weak form of the Ampère equation, i.e. [START_REF] Dular | A posteriori error estimation of finite element solutions via the direct use of higher order hierarchal test functions[END_REF],

(h,curl a ') ! "( j,a ') ! +< n# h,a ' > $ = 0 , ! a ' " F 1 (#) , (1)
where F 1 (Ω) is a curl-conform function space defined on Ω; ( • , • ) Ω and < • , • > Γ denote a volume integral in Ω and a surface integral on Γ, respectively, of the product of their field arguments; n is the unit normal on Γ exterior to Ω. Magnetic field h and electric current density j are related to magnetic flux density b and electric field e (in Ω c ⊂ Ω, with Ω = Ω c ∪ Ω c C ), respectively, through the material relations h = µ -1 b and j = j s + σ e, with µ the magnetic permeability, σ the electric conductivity and j s the source current density in stranded inductors in Ω s ⊂ Ω c C . With the magnetic vector potential a defined via b = curl a and e = -∂ t a -grad v (∂ t is j ω with ω the angular frequency in the frequency domain), weak formulation [START_REF] Babuska | A posteriori error estimates for the finite element method[END_REF] gives the a-v formulation (or A-ϕ formulation with notations of [START_REF] Creusé | Residual-based a posteriori estimators for the A/Phi magnetodynamic harmonic formulation of the Maxwell system[END_REF][START_REF] Dular | A posteriori error estimation of finite element solutions via the direct use of higher order hierarchal test functions[END_REF]), i.e.

(µ !1 curl a,curl a ') ! +("# t a,a ') ! c +(" grad v,a ') ! c !( j s ,a ') " s +< n# h,a ' > $ = 0 , ! a ' " F 1 (#) , (2)
with F 1 (Ω) gauged in Ω c C , and containing the basis functions for a as well as for the test function a' (at the discrete level, this space is defined by edge FEs; the gauge is based on the tree-co-tree technique); note that F 1 (Ω) can be non-gauged if a proper resolution scheme is used [START_REF] Creusé | Residual-based a posteriori estimators for the A/Phi magnetodynamic harmonic formulation of the Maxwell system[END_REF]. Once a solution is obtained in a particular discrete function space with (2), it can be used a posteriori to estimate the associated error.

B. Residual error estimator (REE)

Generally, the discretization error committed by the FE method can be resumed by the error on the energy. In the case of the a-v formulation, the local error on a and v can be written e a = aa h and e v = v -v h , with a and ν the exact solutions of the problem and a h and v h the numerical ones. From these local errors, the global error norm can be written, involving error terms on magnetic energy and Joule losses, i.e.

! 2 = µ -1 " # curl e a 2 + $ % " c # j%e a + grad e v 2 . (3)
Carrying out an ad-hoc Helmholtz decomposition of (1) and using suitable interpolation operators and Cauchy-Schwarz inequality, the efficiency and the reliability of REE can be obtained [START_REF] Creusé | Residual-based a posteriori estimators for the A/Phi magnetodynamic harmonic formulation of the Maxwell system[END_REF]. For the a-v formulation, the REE η ω on FE ω ⊂ Ω quantifies all the terms weakly verified in (2) via

! " 2 = ! ";1 2 + ! ";2 2 + ! ";3 2 + (! #;1 2 + ! #;2 2 ) #$%" & (4)
with

! ";1 = h " # h j s $ curl(µ $1 curl a) $ %( j"a + grad v) " , (5a)
! ";2 = h " div(#( j"a + grad v)) " , (5b)
! ";3 = h " j s # $ h j s " , (5c)
! ";1 = h " 1/2 [n #µ $1 curl a] " "' , (5d)
! ";2 = h " 1/2 [#( j$a + grad v) % n] " "' . ( 5e
)
Operator π h represents the projection operator from H(div, Ω) in a discrete approximation space. Coefficients h ω and h γ stand for the diameter of the smallest sphere containing the element ω and the diameter of the smallest circle containing the facet γ, respectively. All these terms must be evaluated on each FE. Terms (5a-b) quantify the error of the equation to solve. Term (5c) corresponds to the discretization error of the source term. Finally, terms (5d-e) evaluate the jumps (discontinuities) of the tangential component of h and the normal component of j through the element facets.

C. Hierarchical error estimator (HEE)

Once weak form (2) has been solved, it is evaluated with the corresponding solution with HO basis functions (BFs) a p ' ∈ F p 1 (Ω) as TFs. This gives a left hand side of ( 2) not equal to zero anymore but defining a residue R p that can serve as an error estimator, called HEE [START_REF] Dular | A posteriori error estimation of finite element solutions via the direct use of higher order hierarchal test functions[END_REF], i.e.

R p = (µ !1 curl a,curl a p ') " +(#$ t a,a p ') " c +(# grad v,a p ') " c !( j s ,a p ') " s +< n# h,a p ' > $ , a p ' ! F p 1 (") , (6)
Hierarchical HO-BFs can be associated with edges, facets and volumes of finite elements (FEs) [START_REF] Webb | Hierarchal scalar and vector tetrahedra[END_REF]. For an edge or a facet HO BF a p ', respectively in 2-D and 3-D, shared by two FEs ω a and ω b , with interface γ ab , residue (6) has a contribution only from these two FEs. Restricting the form [START_REF] Webb | Hierarchal scalar and vector tetrahedra[END_REF] to ω a and ω b separately, with a p ' equal to zero on all the boundaries except γ ab , integrating their first volume integral term by parts, and adding the resulting forms, one gets

R p = (curl h ! j,a p ') " a +" b !< n ab # (h | " a ! h | " b ,a p ') > $ ab , (7)
with n ab = n| ∂ω a = -n| ∂ω b the unit normal to γ ab exterior to ω a . Residue R p , calculated via ( 6), is therefore a union of two errors, pointed out by ( 7): an error on the equation to solve, via the volume integral term on ω a + ω b , and an error on the tangential continuity of h, via the surface integral term on γ ab . For an HO-BF a p ' associated with the volume of a FE ω a , thus with a zero trace on ∂ω a , the residue given by ( 6) is reduced to

R p = (curl h ! j,a p ') " a . ( 8
)
Using TFs a' = grad v' in (2) (which is implicitly the case with edge FEs for a' , that contains the gradient of nodal FEs), weak form (2) becomes the weak form of div j = 0, i.e.

(!" t a,grad v ') # c +(! grad v,grad v ') # c $< n % j,v ' > & s =0 , (9)
with Γ s some possible cross sections of conductors in Ω c . A residue R v,p can be also defined for (9) when using HO TFs v p ' for v', being also equal to

R v, p = (div j,v p ') ! a +! b "< n ab # ( j | ! a " j | ! b ,v p ') > $ ab . ( 10
)
Residue R v,p points out an error on div j = 0 on ω a + ω b and an error on the normal continuity of j on γ ab . A useful additional step can consist in using R p as a source (right hand side) for a local FE problem (limited to the FE support of each TF), calculating the HO correction a p to be given to solution a for satisfying the HO weak formulation [START_REF] Dular | A posteriori error estimation of finite element solutions via the direct use of higher order hierarchal test functions[END_REF], i.e.

(µ !1 curl a p ,curl a p ') " +(#$ t a p ,a p ') " c +(# grad v p ,a p ') " c =! R p . ( 11
)
This allows an error estimation directly relative to the scale of the solution, whereas R p is only relative to the others residues for the same mesh (as for the REE; the direct use of R p can define another kind of REE).

III. COMPARISON OF ERROR ESTIMATORS

REE and HEE show some equivalences and differences. Equivalences can be found in all the volume and surface integral contributions involved in the estimators, where both weakly defined equations and interface conditions are evaluated in a way to express their related contributions to the errors. A difference appears regarding the weighting of the contributions, with coefficients defined as smallest diameters for the REE while such coefficients are implicitly defined via integrations in each FE with the HEE. This gives differences, in particular in anisotropic meshes.

Local contributions to the errors are related to each single FE with the REE and with groups of FEs (sharing edges or facets) or single FE (for volume TFs) with the HEE. In both estimators, the contributions are only relative to other similar contributions (which can be useful for mesh adaption), and are thus not expressed at the scale of the solution. However, the HEE can be extended to give a solution-relative error by calculating HO corrections.

Fig. 1 .

 1 Fig. 1. Example of distribution of local error estimation in an eddy current problem. The highest error is located in the skin depth with a too coarse mesh.