Guillaume Losfeld

Vincent Escande

Tanguy Jaffré

Laurent L'huillier

Claude Grison
email: claude.grison@cefe-cnrs.fr

The chemical exploitation of nickel phytoextraction: An environmental, ecologic and economic opportunity for New Caledonia

Keywords:

published or not. The documents may come

Introduction

New Caledonia is estimated to account for 20-25% of the world nickel resources. Industrial activities around nickel account for 10% of the island's gross domestic product and could increase to 30-40% in the upcoming years [START_REF] L'huillier | Mines et environnement en Nouvelle-Calédonie : les milieux sur substrats ultramafiques et leur restauration -Editions IAC[END_REF]. New Caledonian mines are opencast and before exploitation starts, vegetation, topsoil and upper soil horizons need to be removed. As high-grade ores such as garnierite are becoming scarce, new low-grade lateritic ores are being used with larger excavations. This inevitably leads to increased waste production in ore processing [START_REF] L'huillier | Mines et environnement en Nouvelle-Calédonie : les milieux sur substrats ultramafiques et leur restauration -Editions IAC[END_REF]. Thus the overall impact of mining on ecosystems is serious, with perturbations on water flows, soil erosion, contaminated sediment transport and bioaccumulation issues [START_REF] Fernandez | A combined modelling and geochemical study of the fate of terrigenous inputs from mixed natural and mining sources in a coral reef lagoon (New Caledonia)[END_REF][START_REF] Hedouin | Nickel bioaccumulation in bivalves from the New Caledonia lagoon: seawater and food exposure[END_REF][START_REF] L'huillier | Mines et environnement en Nouvelle-Calédonie : les milieux sur substrats ultramafiques et leur restauration -Editions IAC[END_REF]. Technical feasibility and costs obviously limit soil remediation initiatives: in the case of metal pollution, the conventional method is excavation followed by burial at a waste site for an estimated cost of $400,000 per treated hectare [START_REF] Raskin | Phytoremediation of toxic metals: using plants to clean up the environment[END_REF]. Longterm solutions with lower operating-costs are increasingly considered and for metals, phytoremediation appears as a sustainable technique to reclaim land. It is aesthetically more pleasing and boasted better public acceptance where used [START_REF] Singh | Phytoremediation: an overview of metallic ion decontamination from soil[END_REF]. Phytoremediation has a high potential in New Caledonia: the island is a metallophyte flora hotspot with 2153 species (83% endemics) identified on ultramafic soils [START_REF] L'huillier | Mines et environnement en Nouvelle-Calédonie : les milieux sur substrats ultramafiques et leur restauration -Editions IAC[END_REF]. 27 past experiences in phytoremediation have been reported: these programs mainly considered species selection and agronomic techniques for the development of a sustainable plant cover on mine sites [START_REF] L'huillier | Mines et environnement en Nouvelle-Calédonie : les milieux sur substrats ultramafiques et leur restauration -Editions IAC[END_REF]. Nickel phytoextraction was not considered, because a considerable limit of phytoextraction remained: without credible outlets to dispose of metal-rich biomass, phytoextraction only provides little improvement [START_REF] Sas-Nowosielska | Phytoextraction crop disposal -an unsolved problem[END_REF]. The use of nickel-rich biomass to produce catalysts used in organic syntheses could bring about change and spur the development of phytoextraction in New Caledonia.

Chemical exploitation of metal-rich biomass

Nickel recovery from biomass

The first step of this project was to develop an efficient method to recover metals from metal-rich biomass. Our main concern was to avoid the use of organic solvents or separation agents as they are hazardous and costly substances. Thus the selection of conditions guided us towards a simple, low-cost and efficient process in accord with the principles of green chemistry. It is also compliant with direct industrial applications [START_REF] Grison | Use of metal-accumulating plants for the preparation of catalysts that can be used in chemical reactions[END_REF].

Leaves of Psychotria douarrei and Geissois pruinosa were harvested in the South province of New Caledonia. Biomass is dried at 60 °C (until a stable weight is achieved) and crushed. A thermal treatment at 500 °C (5 h) is then applied in order to destroy organic matter. The ashes obtained are treated with hydrochloric acid to finish the destruction of the remaining organic compounds and to convert metallic cations into metal chlorides. In a typical experiment, the reaction mixture was stirred for 2 h at 60 °C, and then filtered on Celite. The resulting solutions, composed of different metal chlorides, were then concentrated and dry residues could be stored in a stove at 90 °C for several weeks without undergoing hydrolysis. Table 1 summarizes the main steps of the process and indicates mass variations.

20 mL 1 M HCl per gram ashes allowed total recovery of the metals present in the plant ashes. Inductively Coupled Plasma-Mass Spectrometry (ICPS-MS) analyses were used to determine the composition of the plant extracts obtained. Table 2 exemplifies the composition of the different plant extracts obtained.

Use of metal hyperaccumulator biomass in organic chemistry

In a chemical reaction, a catalyst should be added in small quantities compared to other reagents in order to increase the reaction kinetics. As a general rule, 2-5% metals are enough to convert 1 mol of substrate. Using metal hyperaccumulators' biomass we also expected to benefit from the natural polymetallic composition of the plant extracts obtained [START_REF] Corma | Lewis acids: from conventional homogeneous to green homogeneous and heterogeneous catalysis[END_REF]. The use of catalysis was established as the 9th of the twelve principles of green chemistry and it is sometimes referred to catalysis as the ''foundational pillar'' of green chemistry [START_REF] Anastas | Catalysis as a foundational pillar of green chemistry[END_REF]. Heterogeneous catalysis presents itself as the latest improvement as it allows catalyst and reagent recycling: active species (metal hyperaccumulators extracts here) are dispersed on a support (montmorillonite K10 here).

The acetylation of anisole (Fig. 1) in heterogeneous catalysis is a witness of the recent evolution of chemistry towards new, more environmentally friendly processes: the reaction has been extensively studied and is a very interesting model to test the efficiency of 'green' Lewis acid supported catalysts. This choice was also motivated by the industrial interest of methoxyacetophenone. Our results using nickel hyperaccumulators extracts supported on montmorillonite K10 are summarised in Table 3. The increased activity of NiCl 2 in heterogeneous catalysis is very interesting considering the advantages of heterogeneous catalysis over homogeneous catalysis (Entry 2-3).

The following features should be noted: (i) The catalysts derived from P. douarrei and G. pruinosa were thermally stable and easy to use. (ii) The products could be isolated without work-up. (iii) A simple filtration and thermal activation allowed the green catalyst reuse. (iv) The green catalysts kept a high activity after four recycling (Table 4). (v) The green catalyst activity was much higher than commercial NiCl 2 . (vi) The green catalyst efficiency was not affected by the formed acetyl moiety, which is a great advantage for Friedel-Crafts acylation. Smaller amounts of green catalysts were required, thus an improved environmental quotient (EQ) was obtained. It also appeared that in heterogeneous catalysis the para/ortho ratio increased dramatically (up to 100% para isomer). This is due to a rearrangement of the product towards the more thermodynamically stable para compound: it is a positive effect as it avoids purification issues. The best results were obtained using a P. douarrei extract after solvent-less purification: the Dowex TM M4195 (The Dow Chemical Company) ion exchange resin was used to concentrate nickel: 100% yield was achieved in 6 h at 70 °C.

This shows that our process is efficient and relevant in today's chemistry context: it performs better than other cutting-edge heterogeneous catalytic systems [START_REF] Spagnol | Selective etherification of phenylketone cpds. in 4-position|part. suited for the prepn. of 4-methoxyacetophenone from anisole[END_REF][START_REF] Gupta | Silica supported zinc chloride catalyzed acetylation of amines, alcohols and phenols[END_REF]. Moreover, the source of metals should also prove less costly than conventional metals obtained from mining.

Comparison with existing outlets for metal-rich biomass

To date, two outlets for metal-rich biomass have been thoroughly investigated: 'phytomining' (Brooks et al., 1998), an alternative metal production process developed at the United-Stated Department of Agriculture (USDA) for nickel and bio-energy production [START_REF] Van Ginneken | Phytoremediation for heavy metal-contaminated soils combined with bioenergy production[END_REF]. Other minor outlets exist, such as compost production to supplement metal deficient soils [START_REF] Tang | Effects of copper-enriched composts applied to copper-deficient soil on the yield and copper and zinc uptake of wheat[END_REF] or metal nanoparticle production [START_REF] Qu | Zinc accumulation and synthesis of ZnO nanoparticles using physalis alkekengi L[END_REF].

Bioenergy production

Although comprehensive economic studies showed the potential of bio-energy, a large quantity of biomass is required which leads to restrictions on species used: willow in short-rotation

Table 1

Mass obtained from 100 g dry leaves in each step of the process (g ± standard error).

Dry leaves

Thermal treatment Acid treatment P. douarrei 100 11.2 (±0.4) 15 g (±0.9) G. pruinosa 100 14.5 (±0.4) 9.8 g (±0.5) coppice, rapeseed, maize and wheat have been considered [START_REF] Van Ginneken | Phytoremediation for heavy metal-contaminated soils combined with bioenergy production[END_REF]. Although they do not specifically target any metal, these species were selected for their high annual biomass production and their general tolerance to trace metals. However plant ashes with a low metal content are to be considered as problematic toxic wastes. Moreover such species are not available in New Caledonia: species found on ultramafic soils (including metal hyperaccumulators) are rather slow growing [START_REF] L'huillier | Mines et environnement en Nouvelle-Calédonie : les milieux sur substrats ultramafiques et leur restauration -Editions IAC[END_REF] and are not suitable for bio-energy production.

Nickel phytomining

In phytomining, nickel hyperaccumulators are considered as a 'bio-ore' and the commercial value of the metal produced could balance the costs involved in the development of phytoextraction. It targets highly contaminated soils and low-grade ores where conventional mining techniques would be too costly. As leaves need to be burned to ashes in the nickel recovery process, a possible production of energy is also possible: expected returns in the best case including nickel and energy production are estimated to $1311 per hectare (Brooks et al., 1998). Previous experiences showed that yields could be good enough to foresee commercial applications.

The alternative we propose, which is to use the metal-rich biomass to produce Lewis acid catalysts is conceptually very different from phytomining and could overcome some limits of this technology. It specifically tackles metal contamination regardless of the metal considered and it is also compliant with energy recovery. The end use of metals recovered, which is organic chemistry, allows broader perspectives. Possible valorisation of New Caledonian nickel hyperaccumulators P. douarrei and G. pruinosa is considered but zinc hyperaccumulators could also prove very useful. Metallic species are directly obtained in the required state to produce Lewis acid catalysts: herein, nickel chloride is obtained in a simple three steps process from biomass. Compared to currently-used metal production processes followed by metal-chlorides production, the number of steps involved is lesser, and detailed lifecycle analyses should prove the process we describe less costly. Moreover, as exemplified by the acetylation of anisole, purification is not mandatory: although roughly purified P. douarrei extracts proved more efficient than untreated plant extracts, a mixture of metallic species still showed increased activity compared to pure commercial nickel chloride.

Conclusion

Phytoextraction is currently raising much interest as a sustainable method to remediate metal pollution. However, to be fully applicable this innovative technique needs strong incentives including credible outlets for metal-rich biomass. It appears that existing outlets could not give phytoextraction sufficient momentum to lead to a general use of the technology. The alternative we propose, using metal hyperaccumulator plants species as a source of Lewis acid catalysts for organic chemistry could bring about change: it is more adapted to the state of metallic species present in plant tissues; and it could bring extra added value to the metal-rich biomass produced. In terms of organic synthesis, perspectives are various and include hydrochlorinations of alcohols, synthesis of aromatic heterocycle scaffolds, chiral cyclic structures, protections of alcohols and amines, aldolisations, metal catalysed rearrangements and other key syntheses of industrial and fine chemistry.

Table 2

 2 Composition (mass percentage) of the plant extracts obtained after treatment.

		Mg	Al	Ca	Mn	Fe	Ni
	P. douarrei	1.6	0.3	6.0	0.1	0.4	11.5
	Purified P. douarrei	0.4	0.1	2.6	0.0	0.1	22.3
	G. pruinosa	7.6	0.5	12.0	0.2	0.8	7.4
				Catalyst
	Fig. 1. Acetylation of anisole using acetic anhydride.

Acknowledgements

The authors gratefully acknowledge financial support from the Centre National de la Recherche Scientifique (CNRS), the PIR-CNRS-CEMAGREF INGECOTECH and the ANR (11ECOT01101) and two Grants from the Ecole Polytechnique, Palaiseau and ADEME. They also thank Bruno Buatois for the CEFE platform of chemical analysis, Eddy Petit (IEM/UM2 UMR 5635) for the NMR analysis and Rémi Freydier (Hydrosciences laboratory/AETE platform/UM2) for the ICP-MS analysis.


Frédéric Lebon
email: lebon@lma.cnrs-mrs.fr

Michel Raous
email: raous@lma.cnrs-mrs.fr

Iulian Rosu
email: rosu@lma.cnrs-mrs.fr

Multigrid Methods for Unilateral Contact Problems with Friction

Keywords: Unilateral contact, friction, multigrid

In this paper, multigrid methods are tested on unilateral problems with friction. An optimal strategy is presented and efficiency of the solver is discussed on several examples.

Introduction

Multigrid methods have been widely used in fluids mechanics when large numbers of degrees of freedom are involved. Usually the geometries are sufficiently simple to enable the generation of multiple overlapped meshes in an easy way (essentially in the context of finite difference methods). In nonlinear structure mechanics, the computational costs increase because of the treatment of nonlinearities and finite elements methods are dominant because of the complexity of the geometries. The present work investigates the ability of multigrid methods to reduce the computational times and analyzes the specific problems of formulation and implementation related to the treatment of nonlinearities in the context of finite element methods. This work is conducted on contact problems involving unilateral contact and friction between an elastic body and a rigid obstacle. The nonlinearities are stiff because the contact behavior laws are nonsmooth (the nonpenetration is characterized by the nonregularized Signorini conditions) and nondifferentiable because of the use of the nonregularized Coulomb law.

The Frictional Contact Problem

The Contact Model

Details of the formulation are to be found in [START_REF] Lebon | Multibody contact problem including friction in structure assembling[END_REF][START_REF] Raous | Quasistatic signorini problem with coulomb friction and coupling to adhesion[END_REF][START_REF] Raous | Numerical methods for frictional contact problems and applications[END_REF]. Considering 2D problems, let be an open bounded set of lR 2 , which is the interior of an elastic obstacle, with a sufficiently smooth boundary = F ∪ D ∪ C , in contact with a rigid body. The boundary D (respectively F ) is the part of on which the displacements (respectively the forces) are prescribed. Part C is that on which the contact is possible with the rigid obstacle. Let σ be the stress tensor, K the elasticity tensor, ε the strain tensor and u the displacement.

• Unilateral contact conditions

On the contact boundary C , the displacement u and the unknown contact force F are written in local coordinates (n,t) where n denotes the outward normal unit vector to C :

σ.n = F on C (1) u = u N n + u T ( 2
)
F = F N n + F T (3)
On C , the unilateral conditions are written as follows where d denotes the initial gap:

u N ≤ d (4)
F N ≤ 0 ( 5 )

(u N -d).F N = 0 ( 6 )
This is a Signorini problem. Neither regularization nor penalty are used.

• The Coulomb friction law It is written as follows where µ denotes the friction coefficient:

F T ≤ µ|F N | with if F T < µ|F N | then uT = 0 ( 7 ) if F T = µ|F N | then uT is colinear
and opposite to F T [START_REF] Duvaut | Inequalities in Mechanics and Physics[END_REF] As presented in [START_REF] Cocu | Formulation and approximation of quasi-static frictional contact[END_REF][START_REF] Raous | Quasistatic signorini problem with coulomb friction and coupling to adhesion[END_REF][START_REF] Raous | Numerical methods for frictional contact problems and applications[END_REF], by introducing an incremental formulation, it was established that the problem to be solved at each step, is subsequently very similar to the static one. In that case, the unknowns are the displacement increments and there are extra terms both in the functional which has to be minimized and in the convex characterizing the constraints. In this paper, we consider only the static problem using the following displacement formulation for the friction (( 9) and ( 10)). Formulation and numerical algorithms can easily be extended to quasistatic evolutive problems [START_REF] Cocu | Formulation and approximation of quasi-static frictional contact[END_REF][START_REF] Raous | Quasistatic signorini problem with coulomb friction and coupling to adhesion[END_REF]. As for the unilateral contact, no regularization (compliance) is used.

F T ≤ µ|F N | with if F T < µ|F N | then u T = 0 ( 9 ) if F T = µ|F N | then u T is colinear and opposite to F T (10)

Variational Formulation

The variational formulation of the Coulomb problem can be written as an implicit variational inequality [START_REF] Duvaut | Inequalities in Mechanics and Physics[END_REF]. By introducing a fixed point problem P 1 on the sliding limit g [START_REF] Panagiotopoulos | Inequality Problems in Mechanics and Application[END_REF][START_REF] Raous | Quasistatic signorini problem with coulomb friction and coupling to adhesion[END_REF], the problem can be solved as a sequence of classical variational inequalities including nevertheless a nondifferentiable term. Because of the symmetry of the elasticity mapping, these problems can be solved as minimization problems under constraints (problem P 2 ).

Problem P 1 : Find the sliding limit function G, fixed point of the application g → µ|F N (u g )| where u g is solution to the problem P 2 depending on g.

Problem P 2 : Let φ 1 be the load applied in , φ 2 be the load applied to part F of the boundary, and g be given,

find u g ∈ lK = {v, v = 0 on D and v N ≤ d on C } such that: J (u g ) ≤ J (v) ∀ v ∈ lK (11)
with

J (v) = 1 2 a(v, v) -(f, v) + j (v) ( 12
)
where: a(u, v) = grad s u K grad s v dx (13)

j (v) = C g v T dl (14) (f, v) = φ 1 v dx + F φ 2 v dl (15)

Preconditioner for Mixed Formulation

In this section, we adopt the mixed formulation proposed in [START_REF] Alart | A mixed formulation for frictional contact problems prone to newton like solution methods[END_REF]. The problem is split into two parts involving the pair x = (u, λ), where u denotes the displacement vector and λ the contact nodal forces, say a differentiable part G and a nondifferentiable one

F : G(x) + F (x) = 0 ( 16
)
This last problem is solved using the Generalized Newton Method (GNM):

x i+1 = x 1 -(K i + J i ) -1 (G(x i ) + F (x i )) ( 17
)
where K i = ∂G(x i ) and J i ∈ ∂F (x i ). The tangent matrix A = K i + J i is nonsymmetric (with zeros on the diagonal) and the previous linear system can be solved by the Conjugate Gradient Squared Method (CGS) associated with the coarse/fine preconditioner proposed in [START_REF] Alart | Solution of frictional contact problem using ilu and coarse/fine preconditioners[END_REF]. The basic idea of the preconditioner is to assume that the tangent matrix A has been obtained through the refinement of a coarser mesh. The matrix A is split into four parts (index c refers to coarse level and f to the fine one):

A = A ff A f c A cf A cc ( 18
)
A block factorization of matrix A gives:

A = A ff 0 A cf A cc -A cf A -1 ff A f c I A -1 ff A f c 0 I ( 19
)
We postulate the following preconditioner:

C = A ff 0 A cf A * cc I A -1 ff A f c 0 I ( 20
)
where A * cc is the tangent matrix on the coarse grid. The coarse/fine preconditioner requires to solve two sub-systems with the matrix A ff . It is possible to define an approximated LDU factorization of this matrix associated to each contact status [START_REF] Alart | Solution of frictional contact problem using ilu and coarse/fine preconditioners[END_REF]. This preconditioner is very efficient because it is not necessary to find an optimal parameter.

Formulation and Implementation of the Multigrid Method for Contact Problems

Classically, in the linear case, the efficiency of multigrid methods can be explained in terms of the spatial frequencies of the error between the solution and the iterates when this error is written in the basis of the eigenvectors of the discretized mapping. The smoothers (Gauss-Seidel method for example) are known to quickly reduce the high frequencies of the error, but a large number of iterations is required to reduce the low components. Multigrid methods increase the convergence on the low frequencies of the error. They correct the solution obtained after only a few smoothings on the fine grid (where the solution is needed) by a correcting error computed on a coarse grid with the current defect of equilibrium. The sequence is then composed of a few smoothings of the solution on the fine grid (large system) and of complete resolution to determine a correction error on a coarse grid (small system). Appropriate interpolation and restriction have to be carried out between the grids.

With linear problems, the problem to be solved on the coarse grid is similar to the initial one defined on the fine grid. With nonlinear problems, one of the difficulties is to write an appropriate problem on the coarse grid, which is defined on the error and associated with the defect of equilibrium computed on the fine grid and carried onto the coarse one. Here, several alternatives arise. Section 4.2.1 is devoted to the problem written on the error. In Section 4.2.2 different methods of carrying the nonlinearity information needed onto the coarse grid from the fine one are proposed and tested.

Basic Solvers

The Discrete Problem

Finite element discretization of problems P 1 and P 2 gives the following problems P 3 and P 4 written in lR 2N h , where N h is the number of fine mesh nodes. Let lK h ⊂ lR 2N h be the cone defined by:

lK h = K h,1 × K h,2 × . .

. × K h,2N h

Let I h be the set of indexes of the equations relating to the normal components of the contact nodes. Let L h be the set of suffixes of the tangential components of the nodes of C .

K h,i = ] -∞, d ] if i ∈ I h K h,i = lR otherwise
Problem P 3 : find the sliding limit G h fixed point of the application g h → µ|F h,N (u h )| where u h is the solution to the problem P 4 depending on g h .

Problem P 4 : find u h ∈ lK h such that:

J h (u h ) ≤ J h (v h ) ∀ v h ∈ lK h ( 21
)
where:

J h (v h ) = 1 2 v T h A h v h -v T h f h + |v h | T g h ( 22
)
with

|v h | T = (|v h,1 |, |v h,2 |, . . . , |v h,2N h |) ( 23
)
|v h | T g h = i∈L h |v h,i |g h,i (24)
and

-A h = (a h,ij ) is the matrix (2N h × 2N h ) associated to the elasticity operator, -f h = (f h,ij )
is the load vector, g h is the vector associated to the sliding limit g (g h ∈ lR card(L h ) ).

The discrete problem P 4 on the fine grid has been derived from an appropriate form of the approximation problem given in [START_REF] Licht | Remarks on a numerical method for unilateral contact including friction[END_REF].

The problem P 3 is solved using a fixed point algorithm (successive approximations). Starting with g 0 h equal to zero (frictionless case), the solution to problem P 4 gives the displacement u h and the contact force F h,N at each step (F h,N is computed directly from the defect of equilibrium A h u h -f h ). We then define the iterative process as g +1 h = µ|F h,N |.

Gauss-Seidel Algorithm for the Contact

Problem P 4 , where the sliding limit g h is given by the fixed point method, is solved by an internal multigrid process. The first step in the multigrid process is the smoothing step for which a Gauss-Seidel method is often used. Here we consider a Cryer-Christopherson method [START_REF] Cryer | The solution of a quadratic programming problem using systematic overrelaxation[END_REF] adapted to deal with the nondifferential term due to friction [START_REF] Lebon | Multibody contact problem including friction in structure assembling[END_REF][START_REF] Raous | Quasistatic signorini problem with coulomb friction and coupling to adhesion[END_REF][START_REF] Raous | Numerical methods for frictional contact problems and applications[END_REF]. It is a Successive Over-Relaxation method with Projection (SORP):

u k+1 h,i = P K h,i ((1 -ω)u k h,i + ωm) (25)
where

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ m = 1 a h,ii ⎛ ⎝ f h,i - i-1 j =1 a h,ij u k+1 h,j - 2N h j =i+1 a h,ij u k h,j -sgn(m)g h,i ⎞ ⎠ if m ∈ lR * m = 0
e l s e (26) with g h,i = 0 if i ∈ L h . The second alternative of (26) occurs when there is no compatibility between the calculated m and its assumed sign, i.e., we obtain a negative (or positive) value of m when its sign (denoted sgn(m)) was taken to be positive (or negative).

P K h,i is the projection operator on K h,i and ω is the relaxation coefficient

(ω ∈ [1, 2 [).
Because of the lack of theoretical estimates, an optimal value of the relaxation coefficient can be only obtained through dichotomic research and this Over-Relaxation method will be used only when the method acts as a solver on the coarse grid. When the method acts as a smoother on the fine grid, ω is set equal to 1 and the smoother is then the Symmetric Gauss-Seidel method with Projection (SGSP) [START_REF] Grego | Méthodes multiniveaux pour des problèmes de contact unilatéral avec frottement[END_REF].

Mathematical Programming Solver and Complementary Formulation

The frictional contact problem introduced in Section 2 can also be formulated, after discretization, as a complementarity problem [START_REF] Chabrand | Numerical methods for solving unilateral contact problems with friction[END_REF][START_REF] Grego | Méthodes multiniveaux pour des problèmes de contact unilatéral avec frottement[END_REF][START_REF] Klarbring | The treatment of problems in contact mechanics by mathematical programming[END_REF][START_REF] Raous | Quasistatic signorini problem with coulomb friction and coupling to adhesion[END_REF][START_REF] Raous | Numerical methods for frictional contact problems and applications[END_REF]).

Problem P 5 : let f be given, find the displacement u and the contact force F such that:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ A * h u h = f * h + F h for i ∈ I h (normal components) ⎧ ⎪ ⎨ ⎪ ⎩ u h,i ≤ 0 F h,i ≤ 0 u h,i F h,i = 0 for i ∈ L h (tangential components) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ u h,i = λ 2 -λ 1 F h,i = 1 2 (ϕ 1 -ϕ 2 ) λ i ≤ 0 ϕ i ≤ 0 i = 1, 2 λ i ϕ i = 0 ( 27
)
Problem P 5 has been reduced to the contact variables by using a condensation procedure described in [START_REF] Raous | Quasistatic signorini problem with coulomb friction and coupling to adhesion[END_REF]. Matrix A * h is the full sized condensed stiffness matrix which turns out to be nonsymmetric because of the friction, and f * h is the condensed loading force vector. The tangential contact displacement has been separated into positive (forward) sliding λ 1 and negative (backward) sliding λ 2 relatively to the unit vector t. This problem is known as a linear complementary problem, the general form of which is written:

w -Mz = q (28) w ≥ 0 (29) z ≥ 0 ( 30
)
wz = 0 (31)
We solve this problem by using the Lemke method which is a direct pivoting algorithm. Details can be found in [START_REF] Chabrand | Numerical methods for solving unilateral contact problems with friction[END_REF][START_REF] Grego | Méthodes multiniveaux pour des problèmes de contact unilatéral avec frottement[END_REF][START_REF] Klarbring | The treatment of problems in contact mechanics by mathematical programming[END_REF][START_REF] Raous | Quasistatic signorini problem with coulomb friction and coupling to adhesion[END_REF].

The Two Grid Method

The Algorithm

We consider a fine mesh (with N h nodes) and a coarse mesh (with N H nodes). We consider the operator (prolongation) p h H : lR 2N H → lR 2N h and the operator (restriction) p H h : lR 2N h → lR 2N H . Starting with u 0 h ∈ lK h , we define an iterative process:

u k h → u k+1 h ( 32
)
• Smoothing on the fine grid We carry out ν smoothing iterates (two or three iterates) of the SGSP ν procedure (Symmetric Gauss-Seidel with Projection). We denote:

u k+1/2 h = SGSP ν ω (u k h , f h , A h ) (33)
We introduce the variable

e k h = u h -u k+1/2 h
, u h being the (unknown) exact solution of problem P 4 . Now, in order to evaluate the correction, we have to solve the new problem:

Problem P 6 : find e k h (such that e k h + u k+1/2 h
∈ lK h ) such that:

J C h (e k h ) ≤ J C h (w h ) ∀ w h with w h + u k+1/2 h ∈ lK h ( 34
)
where

J C h (w h ) = 1 2 w T h A h w h -w T h (f h -A h u k+1/2 h ) + |w h + u k+1/2 h | T g h ( 35
)
Problem P 6 now has to be set on the coarse grid. Let

r k h = f h -A h u k+1/2 h and r k H = p H h r k h .
• Coarse grid correction Three difficulties occur in treating the nonlinearities of contact problems with friction:

-how to carry the nonpenetration condition from the fine onto the coarse grid, -how to treat the absolute value on the coarse grid, -how to transfer the sliding limit g h from the fine grid onto the coarse one? These points will be discussed in Section 4.2.2. We take U k H a vector chosen in lR 2N H and we set:

lK C H = K C H,1 × K C H,2 × . . . × K C H,2N H with K C H,i = lR if i ∈ I H K C H,i = ] -∞, d -U k H,i ] if i ∈ I H (36)
We define two vectors G H ∈ lR card(L H ) and V H ∈ lR 2N H . As with classical nonlinear multigrid methods, these vectors are chosen to approximate the values on the coarse grid. A new problem (similar to eq. ( 35)) is written on the coarse grid:

Problem P 7 : find e k H ∈ lK C H such that:

J C H (e k H ) ≤ J C H (w H ) ∀ w H ∈ lK C H ( 37
)
where

J C H (w H ) = 1 2 w T H A H w H -w T H r H + |w H + V k H | T G H . How to choose U k H , V k H
and G H will be discussed in Section 4.2.2. Problem P 7 involving the error e k H has to be solved (the loads are replaced by the defect r k H in the second member). With linear problems, it is natural to use a direct method (Cholesky). With contact problems, the SORP method or the Lemke method (linear programming method) will be used. When the Lemke method is used, P 7 must be changed into a complementary problem (P b 7 ) similar to P 5 [START_REF] Grego | Méthodes multiniveaux pour des problèmes de contact unilatéral avec frottement[END_REF][START_REF] Raous | Quasistatic signorini problem with coulomb friction and coupling to adhesion[END_REF]. The solution e H is interpolated onto the fine grid:

e k h = p H h e k H ( 38
)
The previous solution u k+1/2 h is corrected (we must ensure that u k+1 h stays within the convex lK h ):

u k+1 h = P lK h (u k+1/2 h + e k h ) (39)
and we return to Equation (33) until a given degree of precision is achieved.

• Prolongation and restriction operators

The prolongation p h H is a linear interpolation. Consistence of the formulation is insured by choosing p h H = (p H h ) T . With three-node triangles when each element is divided into four triangles, we have

r k H (x) = p H h r k h (x) = r k h (x) + 1 2 y∈C x r k h (y).
C x is the set of nodes connected to x [START_REF] Lebon | Introduction aux méthodes multigrilles en mécanique des structures linéaires et non linéaires[END_REF].

• The initial condition (full multigrid)

The two grid process will be initiated by solving first the problem on the coarse grid (problem P 4 ) and then using the prolongation of this solution u 0 H on the fine grid as the initial condition u 0 h = p h H u 0 H to start the first smoothings. This is known as the full multigrid method. Numerical tests have confirmed that the choice of the full multigrid method turns out to be essential for our problem, as generally observed for nonlinear problems [START_REF] Grego | Méthodes multiniveaux pour des problèmes de contact unilatéral avec frottement[END_REF][START_REF] Lebon | Résolution numérique de problème de frottement de Coulomb[END_REF].

Transfer of the Nonlinearities

• Convex definition

The nonpenetration condition is satisfied by the condition that the solution must be within the convex lK k H in problem P 5 . On the coarse grid, the convex lK k H , associated with the error e k H depends on the solution after the smoothing: it is characterized by the condition (U k H + e k H ) N ≤ d on C . Therefore U k H , which is not needed in linear classical case, has to be evaluated on the coarse grid. There exist various means of evaluating U k H from the computed solution u k+1/2 h on the fine grid. Four of them are considered here:

(a) U k H = q H h u k+1/2 h
where q H h is a restriction. We can use:

-the canonic injection U k H = u k+1/2 h , -the standard operator defined in Section 4.2.1, U k H = 1 2 p H h u k+1/2 h
. Operator p H h has been defined for the forces: it is used here on the displacements, where a coefficient 1/2 has to be introduced, -a more sophisticate operator which involves five nodes of the fine grid. (b) U k H is defined as the largest normal displacement through the node under consideration and its two neighbors:

U k H (x) = u k+1/2 h (z) with z = -argmax {u k+1/2 h (y), y ∈ C x }.
This is a procedure developed by Mandel [START_REF] Mandel | A multilevel iterative method for symmetric, positive definite linear complementarity problems[END_REF]. It is a very good choice because it ensures that u k+1/2 h + e h stays within the convex. This is convenient for theoretical reasons. (c) U k H = u 0 H where u 0 H is the solution initially computed on the coarse grid for the full multigrid initialization defined in Section 4.3. This is derived from the Hackbusch [START_REF] Hackbusch | Multigrid method of the second kind[END_REF] choice for nonlinear problems.

(d) U k H = U k-1 H +e k-1 H (with the initial condition U 1 H = u 0 H )
. This is a variant of (c), where the computed correction e k H is applied at each iteration (Hackbusch with updating). The cases (c) and (d) need the use of the nested iteration or the full multigrid method [START_REF] Hackbusch | Multigrid method of the second kind[END_REF]. The advantage of the full multigrid method is that it provides a good starting solution at each level. The algorithm begins at the coarsest level.

• The undifferentiable term

The problem of the absolute value is very important from the mathematical point of view, and the term j(v) is undifferentiable. The term |e H + V k H | depends on the choice of vector V k H . An approximation of the displacement is needed on the coarse grid. Various choices of V k H are given in [START_REF] Lebon | Résolution numérique de problème de frottement de Coulomb[END_REF].

(a) V k H = q H h u k+1/2 h
where q H h is a restriction. We can use:

-the standard injection V k H = 1 2 p H h u k+1/2 h (tangential components), -the canonic injection. (b) the Hackbusch choice: V k H = u 0 H
where u 0 H is the solution initially computed on the coarse grid for the full multigrid initialization. The use of a full multigrid is necessary in that case. (c) the Hackbusch with updating:

V k H = V k-1 H + e k-1 H (with the initial condition V 1 H = u 0 H ). • The sliding limit
G H is an approximation of the sliding limit on the coarse grid. Two cases are proposed: (a) the restriction of the fine sliding limit at each multigrid iteration, (b) the initial value G 0 H computed with the initial resolution on the coarse grid of the full multigrid method.

Multigrid Method

With very large sized problems, introducing several levels of grids can keep the coarse problem small sized. This is of course important because total resolutions are conducted on the coarse grid. It is easy to generalize the previous algorithm to a MultiGrid Method (MGM) and to implement either V-cycle or W-cycle.

Optimal Strategy: Tests on the Various Alternatives

The multigrid methods have been developed in our own finite element code EU-XENE which is a specific version of GYPTIS dedicated to the contact problems at the LMA [START_REF] Latil | Module gyptis version 1.0: Contact unilatéral avec frottement en mécanique des structures -inéquations variationnelles[END_REF].

The optimal strategy for dealing with the nonlinearities is discussed in the case of a plate which is part of a pressure vessel used in nuclear engineering [START_REF] Lebon | Multibody contact problem including friction in structure assembling[END_REF] (see Figure 1). This is an axisymmetrical problem. The structure is in contact with two rigid solids on an oblique zone (30 • ) and on a horizontal zone. The plate is squeezed between the two rigid obstacles by prescribing a vertical displacement (0.05 mm) of the lower one. Primary and secondary pressures are applied on the two sides of the plate. The numerical tests have been conducted by using five meshes, with 51617 nodes for the finest one and 227 nodes for the coarsest one. Of course, the solution is always computed on the finest one. The meshes are unstructured with local refinements in the contact zone (194 contact nodes on the finest mesh).

Numerical tests have been conducted on this example in order to compare the various alternatives presented in Section 4 for the implementation of multigrid methods for that nonlinear problem. The main results are given and details of the numerical results could be found in [START_REF] Lebon | Multigrid methods for unilateral contact problems with friction[END_REF].

Convex Definition

In this case, the tests are conducted on a frictionless case. The different possible choices for the convex definition turn out to be equivalent. For mathematical reasons, the Mandel process is suitable. This is confirmed by other numerical results obtained on other examples. The Mandel procedure will be used in the computations.

The Undifferentiable Term

The same example is now treated with friction. The convex is defined as chosen in the previous paragraph. We test the various ways of choosing V k H in (P 7 ) as described in Section 4.2.2.

In this case, results have shown that the use of informations from the fine grid can be very unsatisfactory when three grids are implemented. The canonic and standard injection should be avoided. The computations are conducted with 3 pre-smoothings, 2 post-smoothings and by using V-cycle. The "Hackbusch" procedure is always efficient and will be chosen in what follows.

The Sliding Limit

We tested the two alternatives given in Section 4.2.2 for G H in (P 7 ). The computations are conducted with 3 pre-smoothings and 2 post-smoothings. The results have shown that the Coulomb's condition |F T | = µ|F N | is verified with a good accuracy when the "restriction" is used to define the sliding threshold G l H . When the "initial value" is used, the accuracy varies from 5 to 10%. The choice for the coarse sliding limit is therefore the restriction of the fine sliding limit.

Number of Smoothings

Taking the same example, we study the effects of the number of pre-and postsmoothings on the convergence. The results have shown that it is convenient to perform both a few pre-smoothings and a few post-smoothings. We will choose in the following 3 pre-smoothings and 2 post-smoothings which turned out to be convenient choices in the numerical tests.

Type of Multigrid Cycle

Taking the three grid example, we have compared the results using V-cycle and Wcycle on the same example of contact with friction. The W-cycle turns out to be slightly better than the V-cycle when considering the number of resolutions on the coarse grid (number of 3-grid iterates) and the total number of smoothings. In the following W-cycles are used for the case 3-grid or more.

Number of Grids

The number of grids is of course an important parameter. A large number of grids is time and memory size consuming. The tests have shown that in the case of large number of degrees of freedom (very fine mesh) the choice 3 or 4 grids is efficient.

For smaller sized problems 2 grids will be sufficient.

Global Efficiency of the Multigrid Method on a Contact Problem

The efficiency of the method is presented here by comparing the multigrid method to the direct resolution on two examples (other ones could be found in [START_REF] Lebon | Multigrid methods for unilateral contact problems with friction[END_REF]. The example is still that of the axisymmetrical plate in contact with two rigid obstacles (see Figure 1). Results for 26018 and 103234 d.o.f. meshes are given in Table 1. The Lemke solver is used both for the direct resolution and for the resolutions on the coarse grids. The multigrid method divides the total CPU time by 8 to 10. In Figure 2, we have plotted the evolution of the CPU time relative to the number of degrees of freedom. It has to be noted first that control of the efficiency of the multigrid methods for that nonlinear problem was not as simple as for the linear case, because among other features it depends on the contact status. So only general and qualitative comments can be given.

It is interesting to note that in nonlinear structure mechanics such as the present frictional contact analysis, the multigrid method was found to be efficient, even with medium sized discretization instead that for linear problems they are proved to be efficient only for very large sized problems. In the case of a plate with 13009 nodes, the total CPU time was already divided by 8 (Table 1) when the Lemke solver is used.

Figure 2 shows that the rate of convergence does not seem to be strongly modified by the nonlinearity of the problem. In our case no theoretical results exist. In Figure 2, the rate of convergence for the two grid method is quasi linear.

We would like to make some comments on the influence of the computer processor properties on the multigrid efficiency. When we implemented first the multigrid methods, we have tested them on a Dec Alpha processor and also on a CISC processor. Of course, the CPU times are drastically reduced by using the very powerful new processors, but it turned out that multigrid methods were more efficient with the CISC processor than with the recent ones when the multigrid resolution is compared with the direct one. It should be noted that the new generation of processors are highly efficient at computation but still relatively slow in performing addressings and transfers. So simple solvers such as the Gauss-Seidel and SORP ones turn out to be more and more fast because they only do a large number of multiplications and only a small number of transfers. That kind of solvers becomes more and more competitive regarding sophisticated methods such as the Lemke's one and multigrid involving a lot of transfers. So when one refers the direct computation of the solution, it can be observed that efficiency of multigrid decreases when the processor power increases.

Fig. 1 .

 1 Fig. 1. Plate in contact with two rigid bodies.

Fig. 2 .

 2 Fig. 2. CPU time versus the number of degrees of freedom for the vessel pressure plate.

Table 1 .

 1 Direct and multigrid resolution for two meshes of the plate with two contact zones (Lemke solver).

	Number of d.o.f Direct method		2-grid method
		CPU	Iterates Smoothings	CPU
	26018	1 43	6	33	15.7
	103234	28 00	6	3 3	2 16


Julian Erskine
email: julian.erskine@ls2n.fr

Rafael Balderas-Hill
email: rafael.balderashill@ls2n.fr

Isabelle Fantoni
email: isabelle.fantoni@ls2n.fr

Abdelhamid Chriette
email: abdelhamid.chriette@ls2n.fr

Model Predictive Control for Dynamic Quadrotor Bearing Formations

Keywords: Multi-Robot Systems, Formation Control, Model Predictive Control, Quadrotor Swarms, Visual Servoing

Formation control of multi-agent systems deals with groups of robots forming specific spatial geometries. Combined with the advancements of unmanned aerial vehicles (UAVs) in the past decade, formation control may potentially be applied to tasks such as search-and-rescue, surveillance, even collaborative manipulation. A key challenge is the decentralization of formation control, where each agent behaves independently using onboard sensors and computation, improving the scaleability and robustness of the system.

This paper proposes a decentralized controller based on model predictive control (MPC), for the control of formations of quadrotor UAVs defined by inter-agent bearings. The use of MPC allows the controller to account for attitude kinematics, improving upon the results of existing bearing formation control methods based on rigidity and visual servoing approaches, which typically only consider the quadrotor as a single or double integrator. Furthermore the near-optimality of MPC permits a more optimal use of the quadrotors dynamic capabilities for faster maneuvering. Extensive simulations are performed to demonstrate the improved transient formation convergence and fast maneuvering permitted by this controller. Experiments show that it is indeed a real-time feasible solution for bearing formation control.

I. INTRODUCTION

Multi-agent systems (MAS) are of great interest to the robotics community, due to the potential for system redundancy, cost-effectiveness, and augmented capabilities beyond those of individual robots [START_REF] Coppola | A Survey on Swarming With Micro Air Vehicles: Fundamental Challenges and Constraints[END_REF]. A computer able to communicate with all agents can act as a centralized controller, but decentralization is often preferred due to superior scaleability and robustness to communication delays and single point failures. Initial work focused on swarm behaviour and flocking [START_REF] Olfati-Saber | Flocking for Multi-Agent Dynamic Systems: Algorithms and Theory[END_REF], generally following Reynolds' rules [START_REF] Reynolds | Flocks, Herds, and Schools: A Distributed Behavioral Model[END_REF] for which the collective motion is more important than the positioning of individual agents. Recently however there has been significant work in the field of formation control, where precise inter-agent geometries may be optimal for certain tasks such as intruder searches [START_REF] Vencatasamy | Secure a Zone from Intruders with a Group Robots[END_REF], search-and rescue [START_REF] Arnold | Correction to: Search and rescue with autonomous flying robots through behavior-based cooperative intelligence[END_REF], cooperative mapping [START_REF] Kim | Cooperative Robot Exploration Strategy Using an Efficient Backtracking Method for Multiple Robots[END_REF], and even cinematography [START_REF] Mademlis | High-Level Multiple-UAV Cinematography Tools for Covering Outdoor Events[END_REF].

The type of control law for MAS generally depends on the type of information available to the controllers, and on the Fig. 1: Formation flight with three quadrotors constraints of the agents. In situations where the positions of all agents are known, a wide selection of controllers are available [START_REF] Fernando | Formation control and navigation of a quadrotor swarm[END_REF], [START_REF] Kang | Design and Realization of Distributed Adaptive Formation Control Law for Multi-Agent Systems with Moving Leader[END_REF]. The knowledge of all positions in a common frame however is a best-case scenario. Formations of agents that must evolve in unmapped and GPS-denied environments must be able to self-localise with their neighbours using onboard measurements. Often, this takes the form of distance measurements [START_REF] Zelazo | Decentralized rigidity maintenance control with range measurements for multi-robot systems[END_REF], [START_REF] Krick | Stabilisation of infinitesimally rigid formations of multi-robot networks[END_REF] or bearing (often referred to as visual) measurements [START_REF] Montijano | Vision-Based Distributed Formation Control Without an External Positioning System[END_REF]- [START_REF] Zhao | Bearing Rigidity Theory and its Applications for Control and Estimation of Network Systems: Life Beyond Distance Rigidity[END_REF] which return incomplete information about the state of the formation. In this paper we only consider bearing measurements, which are extractable from cheap and lightweight cameras.

Most bearing formation controllers (see [START_REF] Zhao | Bearing Rigidity Theory and its Applications for Control and Estimation of Network Systems: Life Beyond Distance Rigidity[END_REF] for survey) only consider the translation and yaw kinematics of quadrotors, neglecting the full attitude kinematics to which they are coupled (due to the well-known under-actuation of quadrotors), and which become significant during aggressive maneuvering. For agents simulated as single or double integrators this is not an issue, but for real quadrotors this limits the reactivity of the system due to the low gains and filtering which are required to avoid highly dynamic attitudes, and which may result in instability. To solve this, we apply model predictive control (MPC) considering the quadrotors thrust and angular velocity as control inputs to track the desired bearing trajectory of the formation. This allows us to account for the (almost) full dynamics of the quadrotor in the control law, permitting near-optimal convergence and maneuvering of the formation. The major contributions of this paper compared to existing work on MPC for formation control [START_REF] Erunsal | Nonlinear Model Predictive Control for 3D Formation of Multirotor Micro Aerial Vehicles with Relative Sensing in Local Coordinates[END_REF]- [START_REF] Saska | Coordination and navigation of heterogeneous MAV-UGV formations localized by a 'hawk-eye'-like approach under a model predictive control scheme[END_REF] is the use of bearing error in the optimisation cost function, allowing for a weak dependence on the interagent distance. We can thus say our controller is "almostbearing-only" (i.e. functions with only bearing sensors in most situations), while providing a more reactive formation compared to existing bearing-only controllers.

The paper is structured as follows: Section II presents the modelling of quadrotors and bearing formations. In 1 section III, rigidity based control of bearing formations are presented, along with the concept of model predictive control. The application of MPC to bearing formation control is presented in section IV, and is supported by extensive simulations in section V and by real-time experiments in section VI. Finally in section VII we discuss the successes and drawbacks of this novel MAS control strategy.

II. MODELLING

A. Quadrotor Modelling

Each quadrotor i is considered as a rigid body with frame F i = {O i , x i , y i , z i } moving with respect to a common inertial frame F 0 = {O 0 , x 0 , y 0 , z 0 }. The primary interest in decentralized bearing control is to act without any external system of localization. As such, the position p i = ---→ O 0 O i is unknown. The quadrotors are, however, able to measure their velocity v i , acceleration vi , and angular velocity ω i in F i using onboard optical flow sensors, accelerometers, and gyroscopes respectively.

The attitude of the quadrotors are also known. Indeed, an accurate estimation of the roll and pitch (denoted as φ i and θ i ) is necessary for stable flight, and is often recovered from a sensor fusion of the accelerometer and gyroscope. The yaw (denoted as ψ i ) of each agent is more difficult as the magnetic field used for yaw estimations is easily perturbed [START_REF] Svacha | Inertial Velocity and Attitude Estimation for Quadrotors[END_REF], even by the electric current in the UAV itself [START_REF] Sa | System Identification , Estimation and Control for a Cost Effective Open-Source Quadcopter[END_REF]. In the case where a consistent estimate of the absolute yaw is infeasible, a decentralized bearing estimator such as [START_REF] Schiano | A Rigidity-Based De-centralized Bearing Formation Controller for Groups of Quadrotor UAVs[END_REF] can be used for each agent to converge to a common (although arbitrarily rotated) reference frame for the yaw. We choose to use the quaternion q = [q w q x q y q z ] T representation for the attitude control due to it's computationally efficiency.

The state of the i th quadrotor evolves according to:

ṗi = v i (1) vi = f i m i  
2(q w q y + q x q z ) 2(q y q z -q w q x ) 1 -2(q 2 x + q 2 y )

  + g (2) qi = 1 2 0 ω i ⊗ q i (3) ωi = I -1 i τ i -I -1 i (ω i × I i ω i ) (4)
where m i is the mass, f i is the thrust along the z i axis, and τ i is the steering moment of the quadrotor. The inertia tensor I i , as well as ω i and τ i are expressed in F i , while p i , v i , and q i are expressed in F 0 . The gravitational acceleration g = [0 0 -9.81] T ms -2 is also expressed in F 0 . Each quadrotor is controlled by setting the angular velocity of four co-planar propellers. Each propeller p generates a thrust force f p = k t Ω 2 p and a drag torque

τ p = sign(Ω p )k d Ω 2
p , where k t and k d are aerodynamic coefficients and Ω p is the rotational velocity of the propeller. The vector of squared propeller velocities

Ω i = [Ω 2 1 • • • Ω 2 4 ]
T can be directly related to the actuation wrench where Γ depends on the quadrotor geometry and on the aerodynamic coefficients of the propellers [START_REF] Six | Identification of the Propeller Coefficients and Dynamic Parameters of a Hovering Quadrotor from Flight Data[END_REF]. Generally the propeller velocities themselves behave as first order systems but this is neglected along with drag effects and lower level motor and battery dynamics due to difficulties (although not insurmountable [START_REF] Faessler | Differential Flatness of Quadrotor Dynamics Subject to Rotor Drag for Accurate Tracking of High-Speed Trajectories[END_REF]) in accurately identifying them.

w i = [f i τ T i ] T by w i = ΓΩ i (5)

B. Formation Modelling

Typically in formation control, the formation is represented as a graph G(V, E) composed of a set of vertices V and a set of edges E as in Fig. 2a. Each agent A i is represented by a vertex V i ∈ V, and each inter-agent measurement corresponds to an edge

E ij = V i × V j ∈ E. Graphs only represent the flow of information, thus a framework F(G, q) is defined such that q = [q T 1 • • • q T |V| ]
T (known as the embedding) associates a state q i to each V i as shown in Fig. 2b. Because of the under-actuation of quadrotors, the embedding is chosen as the state of the quadrotor projected on the R 3 × S 1 manifold, thus the embeded state of A i is q i = [p T i ψ i ] T . As the quadrotor frame F i evolves in SE(3), a vector x i ∈ R 3 in F i can be reprojected onto the R 3 × S 1 manifold using a ZYX Euler angle convention.

x R 3 ×S 1 i = R y (φ i )R x (θ i )x SE(3) i (6)
We assume that each quadrotor (or agent) is equipped with a camera, capable of measuring the direction (or bearing) of other agents in its own frame of reference (see Fig. 3). The bearing of A j wrt A i expressed in F i can be written as

β ij = R T z (ψ i ) p ij d ij ∈ S 2 (7)
where p ij = p j -p i is the position of A j with respect to

A i expressed in F 0 , d ij = ||p ij ||
is the distance between the two agents, and R z is a rotation matrix about the z 0 -axis. Note that the yaw can be expressed as

ψ = atan2(2q w q z + 2q x q y , -2q 2 y -2q 2 z + 1) (8)
III. CONTROL BACKGROUND To understand the benefits of MPC in the context of formation control, we discuss existing control methods for formations, beginning at the level of the quadrotor and moving up to the higher level formation controller. Fig. 3: Agent A i measuring the bearing β ij . As A i is quasi-static, it has no need to reproject β ij . However A j is performing a high-dynamics maneuver, and thus must reproject any measurements in F j onto its R 3 × S1 manifold.

A. Quadrotor Control Structure

The most common approach to deal with the inherent under-actuation of quadrotors is to have multiple layers of control, with an outer translation controller tracking a highlevel position or velocity trajectory by generating a thrust and reference attitude for the low-level attitude controller. The attitude controller uses the attitude estimation and bodyrate measurements to feed either a kinematic or dynamic controller in order to track the reference attitude by applying actuation moments which, along with the thrust, are then converted to propeller speeds through Eq. ( 5). As a general rule, the high-level planners run at 1-20 Hz, the translation controllers at 20-200 Hz, and the attitude and motor level controllers run at frequencies greater than 200 Hz.

B. Existing Bearing Controllers

The control of bearing formations consists primarily in driving the formation measurements

β = [β T 1 , • • • , β T |V| ] T (where β i (t) = [β ij (t) T • • • β in (t) T ]
T is the stack of bearings observed by A i ) toward to a desired set of values β d . Secondary tasks such as moving the formation may also be imposed. Existing work on quadrotor bearing formations generally uses the first order kinematics of the formation, expressed as

β = B v ψ ( 9
)
where B is the bearing rigidity matrix (similar to the Jacobian in other robotic domains) relating the change in bearings to the velocity

v = [v T 1 • • • v T |V| ] and yaw rates ψ = [ ψ1 • • • ψ|V| ]
of the agents. The quadrotors are then controlled by their body frame velocities and yaw rates such as to minimize the bearing error e β = β dβ. A successful example of this is the rigidity controller from [START_REF] Schiano | A Rigidity-Based De-centralized Bearing Formation Controller for Groups of Quadrotor UAVs[END_REF] which uses the jacobian transpose control law

v ψ = k c B T β d (10)
allowing for a compact, low computation controller which is easily decentralizable and is tuned by a single gain k c . It is important to discuss the previously mentioned steering of the bearing formation. It has been proven in [START_REF] Schiano | A Rigidity-Based De-centralized Bearing Formation Controller for Groups of Quadrotor UAVs[END_REF] that any (rigid) bearing formation in R 3 × S 1 has a five-dimensional matrix in the kernel of the bearing rigidity matrix. This corresponds to the five trivial motions

M = span(ker(B)) = [v F ∈ R 3 ψF ∈ S 1 ṡF ∈ R 1 ] (11)
corresponding to a spatial translation v F of the formation, a rotation ψF of the formation about z 0 , and a change of scale 1 ṡF of the formation. Any trivial motion can thus be expressed in the body-frame velocity of A i as

v i = R T i (v F + ψF (c F -p i ) × z 0 + ṡF (c F -p i )) (12a) ψi = ψF (12b)
where c F is a finite position in space about which the formation rotates and expands, taken in this paper as the centroid of the agents. The fact that the motions in Eq. ( 12) lie in the kernel of the rigidity matrix decouples them from the first order bearing control problem.

Other control laws exist treating agents either as double integrators or as non-holonomic mobile robots, but few formation control works consider the coupling between the attitude and translational dynamics of the quadrotors. In fact, this coupling poses a significant problem for accurate bearing control as the quadrotor is not able to instantly generate the velocity or acceleration set by the formation controller, due to the need to reorient itself. This requires that either the formation control gains are low or that the control signal is heavily filtered, which results in slow convergence to the desired formation and less reactive steering.

C. Model Predictive Control

MPC was initially pioneered to control complex systems with constraints. Recently it has become a popular tool in robotics due to its near optimality (it is truly optimal only for an infinite control horizon [START_REF] Findeisen | An Introduction to Nonlinear Model Predictive Control[END_REF]) and constraint handling capabilities. It has in fact been used in various quadrotor applications trajectory planning and tracking [START_REF] Kamel | Model predictive control for trajectory tracking of unmanned aerial vehicles using robot operating system[END_REF], [START_REF] Neunert | Fast nonlinear Model Predictive Control for unified trajectory optimization and tracking[END_REF], as well as for multi-objective tasks such as tracking a trajectory while keeping and object in view [START_REF] Falanga | PAMPC: Perception-Aware Model Predictive Control for Quadrotors[END_REF].

MPC is performed by predicting the state q of the system over a prediction horizon T p in the future, discretized into N p samples, as a function of the current state q 0 , a desired trajectory q d (t), and the future control inputs u(t). The control input is optimized to minimize an objective function

min u Np i=1 O(q, q d , u, ∆t) s.t. u ≤ u ≤ ū C eq (q, u) = 0 C(q, u) ≤ 0 ( 13
)
where O is the objective function to be minimized, C eq and C are constraint functions which must be respected, and u and ū are the lower and upper bounds of the control signal. Generally once the optimization converges, the first values of the computed control input sequence u(t 0 ) are applied to the system, and the optimization is run again using updated state information.

IV. MPC FOR BEARING FORMATIONS

The first decision is to choose a control variable that may be accurately applied to the quadrotor so that it will evolve as close to the prediction as possible. We choose to use

u i = [f i ω T i ]
T as the quadrotor control, which will be tracked by a low-level control attitude rate controller setting the motor speeds. This is because ω i is measured directly, can be integrated to recover the future attitude, and the attitude rate dynamics of quadrotors are very fast allowing rapid tracking by the low-level rate controller. Motor control is not used because the required rate for stable, vibrationfree flight may be difficult to reliably achieve with onboard real-time optimization, together with potential computation or image processing tasks.

A. Prediction

Given the choice of control signal, we choose the state

q i = [q i , v i , p i ] (14)
for which the prediction q i (t) may be calculated from a control signal u i (t), t ∈ T p and the state q i at t 0 (up to an unknown translation) by integrating Eq. (1-3) with respect to time. Each predicted bearing β ij (t) is then defined by Eq. [START_REF] Mademlis | High-Level Multiple-UAV Cinematography Tools for Covering Outdoor Events[END_REF][START_REF] Fernando | Formation control and navigation of a quadrotor swarm[END_REF]. The position of A j relative to A i expressed in F i given the bearing measurement β ij at time t 0 is predicted as

pj (t) = β ij dij + v F t (15)
where dij is an approximate distance and the v F t term acts as the MPC equivalent of a velocity feed-forward term. Note that this requires that the inter-agent relative yaw is well estimated, which requires communication between connected nodes in the formation graph [START_REF] Schiano | A Rigidity-Based De-centralized Bearing Formation Controller for Groups of Quadrotor UAVs[END_REF]. The distance may be taken as either the desired distance corresponding to a desired formation scale, an estimated distance (using such methods as [START_REF] Schiano | The Dynamic Bearing Observability Matrix Nonlinear Observability and Estimation for Multi-Agent Systems[END_REF], [START_REF] Spica | Active decentralized scale estimation for bearing-based localization[END_REF] for scale estimation), or even a unitary distance as is often done in 2D visual servoing [START_REF] Allibert | Predictive control for constrained image-based visual servoing[END_REF]. As there is no known initial position p i in an inertial frame, the predicted position p i (t) is therefore relative to F i at the time t 0 , projected onto the R 3 × S 1 manifold.

B. Objective Function

The objective function of A i may be a function of the predicted state q i (t), bearings β i (t), and inputs u i (t). As discussed in section III-B, there are two primary tasks of formation control:

1) Achieve the desired inter-agent geometry 2) Steer the formation by it's trivial motions We therefore formulate the objective function as a multiobjective optimization problem of the shape

O i = O βi + wO Mi ( 16
)
where O βi is a function minimizing the bearing error, and O Mi is a function minimizing the trivial motion error, and w is an inter-objective weight.

The bearing error e βi = β d i -β i minimization is achieved using the function

O βi = Np n=1 e T βi (n∆t)Q βi e βi (n∆t) (17)
where Q βi is a positive diagonal matrix.

The second task minimizes the trivial motion error

e Mi = [v d i (t) T ψ d i (t)] T -[v i (t) T ψ i (t)] T (18)
where v d i (t) and ψ d i (t) come from the mapping of a desired trivial motion M d (given by a pilot or a higher level control algorithm for example) to the quadrotor motion by Eq. ( 12). The objective function is thus

O Mi = Np n=1 e T Mi (n∆t)Q Mi e Mi (n∆t) ( 19
)
where Q Mi is a positive diagonal matrix. Because M ∈ ker(B) applies only to the first order model whereas we use a third order model (i.e. the input ω i must be thrice integrated to get the bearing output). An individual minimization of the two objectives may yield conflicting inputs, thus w is used to shape the Pareto frontier of the two tasks.

C. Constraints

The optimization of the objective function may account for constraints that must necessarily be satisfied. There are many possible constraints for this system, but we only briefly discuss them for the sake of completeness.

1) Operational Constraints: These are constraints relating to the quadrotor itself, its actuation limits, operational stability, or similar criteria. The only critical one in our case is the input bounds from Eq. ( 13), which are defined (for our quadrotor) as

f i ≤ f i ≤ fi (20a) ω xy,i ≤ ω xy,i ≤ ωxy,i (20b) ω z,i ≤ ω z,i ≤ ωz,i (20c)
where x and x are respectively the lower and upper limits of x. These constraints are chosen such that there is no motor saturation during an aggressive manual flight respecting these limits. Other constraints could include limitations on the velocity, and on the maximum tilt angle of the quadrotor.

2) Task Constraints: There are no task constraints implemented currently, however we may easily envisage constraints ensuring that all observed quadrotors remain in the field of view of the observing quadrotor as in [START_REF] Schiano | Bearing rigidity maintenance for formations of quadrotor UAVs[END_REF], which require the persistent observation of a target outside the formation [START_REF] Lissandrini | Cooperative optimization of UAVs formation visual tracking[END_REF], or ensuring collision avoidance with environmental obstacles [START_REF] Penin | Vision-Based Reactive Planning for Aggressive Target Tracking while Avoiding Collisions and Occlusions[END_REF].

V. SIMULATIONS

Because of the difficulty in proving the stability of decentralized MPC in this application, this section shows the results of extensive simulations to both demonstrate the stability of the proposed MPC control strategy, as well as to compare its performance relative to other state of the art bearing control methods. We compare our controller with the rigidity-based controller in [START_REF] Schiano | A Rigidity-Based De-centralized Bearing Formation Controller for Groups of Quadrotor UAVs[END_REF], which considers only firstorder quadrotor kinematics, but follows a locally optimal gradient due to the controller of A i considering not only it's observations β ij , but also all observations β ki of itself. In all simulations and experiments, in order to compare results between formations of different sizes, the bearing error magnitude is scaled by 1/ |E|, and thus always has a maximum value of 2.0, regardless of the number of edges. Table I can be referred to for a comparison between bearing magnitude error and the more interpretable angular error.

The quadrotors are modelled in Simulink using Eqns. 1-4 with noise and first order propeller dynamics. The velocity (P), attitude (PD) [START_REF] Brescianini | Nonlinear Quadrocopter Attitude Control[END_REF] and body rate (P) controllers run at 200 Hz, while the outer loop bearing formation controllers run at 30 Hz. We use the well-known Matlab function fmincon to solve the optimization, while the prediction uses a fourth-order Runge-Kutta (RK4) algorithm.

With all agents starting at random positions and yaw angles, and the desired bearings based on another set of random position and yaw angles, the simulation is run for five seconds using both the proposed MPC controller and the rigidity controller. The results of 1110 simulations of the three, four, and five-agent formations (Fig. 4) are shown in Fig. 5a-b, with only formations considered as potentially converging (an average bearing error of e α ≤ 45 • at 5 s) are kept, with the others being assumed to have either crashed (due to controller instability) or fallen into local minimas. It is seen that despite the lack of formal proof, the system shows good convergence properties as only 4% of MPC simulations failed to converge, while 12% of Rigidity controllers failed.

VI. EXPERIMENTS

A. Real-time MPC Implementation

We used the ACADO toolbox [START_REF] Houska | ACADO toolkit -An opensource framework for automatic control and dynamic optimization[END_REF] to generate efficient Ccode for an SQP (sequential quadratic programming) minimization of the objective function with multiple shooting discretization, RK4 integration, and using the QPoases library for solving sparse quadratic programming problems. Since ACADO is pre-compiled, it lacks the runtime flexibility to re-dimension the MPC problem for the variable number of graph edges leaving each vertex (as in Fig. 4b-c). Rather than pre-compile separate code for each agent, we generated code based on a maximum number of observations per agent (we used four) and set the weights of unused measurements to zero at runtime. After testing, the best results were obtained when running the MPC at 50 Hz, performing two steps of the SQP algorithm (about 8 ms per step) each loop and a prediction horizon of 1.5 s discretized into N p = 15 steps. An inter-robot distance estimate is required (see Eq. ( 15)) and we use a constant value of dij = 2.0 m for all measurements.

B. Hardware

The quadrotors are based on the Lynxmotion Crazy2fly 34 cm frame with MT2208 1100 kV brushless motors, 12 A ESCs, 8" × 4" plastic propellers, a 3-cell LiPo battery, and have masses of around 1.0 kg. They are controlled by a Pixhawk running PX4v1.10.1 [START_REF] Meier | PX4: A node-based multithreaded open source robotics framework for deeply embedded platforms[END_REF] which performs sensor fusion and tracks the commands (thrust, body rate, etc...) sent by the formation controller. The PX4 controller was tuned for reactive body rate control with aggressive PD gains, and the integrator removed to reduce phase shift. The formation controller for a given quadrotor runs on its onboard Raspberry Pi 4B computer with 4GB of RAM running Ubuntu 18.04. The ROS Melodic Mavros package is used 4a to handle message exchanges between the computer and the Pixhawk. Each quadrotor receives desired bearing references from the ground station computer over a 5 GHz wireless network. A Qualisys motion capture (MoCap) system with eight cameras is used to send the pose of the quadrotors to the other quadrotors through the ground station. This information is only used for extracting bearing information where needed (using Eq. ( 7), mimicking an onboard monocular camera), and to generate some nullspace commands to assist the pilot.

C. Experimental Procedures and Results

All experiments began with the quadrotors hovering at prespecified positions, and the formation controller takes over at time t = 0 s as seen in Figures 67. As experiments were performed in a small (relative to the formation footprint) 4 × 6 × 4 m flying arena, an autopilot algorithm running at 20 Hz provides the v Fx , v Fz , and ṡF trivial motions in F 0 to keep the formation centered along the short axes or the arena and at a scale of 1 m. The other trivial motions v Fy and ψF were given manually at 50 Hz using a joystick.

Two experiments each are shown for the three (E3.1-2) and four (E4.1-2) agent formations in Fig. 4. We see in figures 6-7 that despite differing pilot inputs, the experiments generally share very similar convergence properties. The difference in E3.1-2 for the initial formation shows that due to slight differences in the initial configuration, the two tests followed different convergence paths and were still successful. We can furthermore see that the maximum steering velocity (2.5 ms -1 in E3.1 and 5.0 ms -1 in E3.2) had some effect on the bearing control, as the peak bearing error for E3.2 at 8 s, 24 s, 49 s and 52 s occur directly after the peak velocity times. This is possibly due to propeller saturation that occurs when attempting to re-orient the quadrotor while maintaining high thrust to reach accelerations of vi > 12 ms -2 .

The second pair of experiments show that similar results may be obtained when scaling up the formation and using variable numbers of edges. Experiments E4.1-2 had a one agent with a single bearing measurement, two agents with two bearing measurements, and a single agent with three bearing measurements. The steering was less aggressive (with peak velocities of 2.5 ms -1 and peak accelerations of 6 ms -2 ) as the formation occupied more volume making it harder to pilot in the arena, and there was more aerodynamic interference between quadrotors. We can see in Fig. 7 the yaw rate of the formation (e.g. 30-36 s in E4.2 among others) visible due to the distinctive sine wave induced by ω zF . A video of the experiments can be found online 2 .

VII. CONCLUSION

This paper showed that model predictive control (MPC) may be applied to bearing formations of quadrotors, improving the transient rate of formation convergence and allowing for high-speed maneuvering of the formation. Extensive simulations have shown that it is able to increase the performance, and possibly the reliability, of bearing-based formations. Future work will include the implementation of MPC formation controllers in constrained environments requiring dynamic maneuvering, and with field of view handling (to which MPC is well suited). Work is underway to integrate real-time visual detection using monocular fish-eye cameras on the quadrotors for the extraction of bearings.

  An embedding of (a) in R 3 × S 1

Fig. 2 :

 2 Fig. 2: The graph and framework representations of a directed bearing formation, where an arrow from A i to A j indicates a measurement β ij

Fig. 4 :

 4 Fig. 4: Test Formation Graphs

  (a) Probability density histograms of the bearing error magnitude of converged simulations at intervals of 1 s (b) A continuous-time box and whisker plot of the bearing error in successful simulations

Fig. 5 :

 5 Fig. 5: The results extensive simulations comparing the MPC (red) and rigidity (blue) controllers.

  The top plot shows the evolution of the normalized bearing error for E3.1 (blue) and E3.2 (orange). The middle and bottom plots shows the desired formation velocity v d F,y (dark green line) and the measured velocities vi,y (light green lines) expressed in F0 for E3.1 and E3.2. respectively. The bearing error peaks correspond to new desired formations (a)-(c).

Fig. 6 :

 6 Fig. 6: Two experimental results with the graph of Fig. 4a

  The top plot shows the evolution of the normalized bearing error for E4.1 (blue) and E4.2 (orange). The middle and bottom plots shows v d F,y (dark green line) and the measured vi,y (light green lines) expressed in F0 for E4.1 and E4.2. respectively. The bearing error peaks correspond to new desired formations (a)-(c).

Fig. 7 :

 7 Fig. 7: Two experimental results with the graph of Fig. 4b

TABLE I :

 I Relationship between bearing error norm ||e β || and the equivalent angular error e α between the measured and desired values of a single bearing measurement

	eα (deg) 180 90	45	30	15	10	5	1
	||e β ||	2.0	1.41 0.77 0.52 0.26 0.17 0.09 0.02

The scale s F is defined as the mean distance between each quadrotor, and the centroid of the formation c F

https://youtu.be/KCpXjPI6eyU

J. Erskine is the corresponding author. This work was supported by the French region Pays de la Loire through the RFI AtlanSTIC2020 project RAPID and by the Natural Sciences and Engineering Research Council of Canada