A Hindawi
email: hindawi@unice.fr

L Rifford
email: rifford@unice.fr

J.-B Pomet

Mass transportation with LQ cost functions

We study the optimal transport problem in the Euclidean space where the cost function is given by the value function associated with a Linear Quadratic minimization problem. Under appropriate assumptions, we generalize Brenier's Theorem proving existence and uniqueness of an optimal transport map. In the controllable case, we show that the optimal transport map has to be the gradient of a convex function up to a linear change of coordinates. We give regularity results and also investigate the non-controllable case.

Introduction

The optimal transport problem can be stated as follows: given two probability measures µ 0 and µ 1 , defined on measurable spaces X and Y respectively and a cost function c : X × Y -→ IR ∪ {+∞} find a measurable map T : X -→ Y which pushes forward µ 0 to µ 1 , that is When the transport condition T ♯ µ 0 = µ 1 is satisfied, we say that T is a transport map, and if T minimizes also the cost, that is if cost c (T ) = min

T ♯ µ 0 = µ 1 (i.
S ♯ µ 0 =µ 1 cost c (S)
then we call it an optimal transport map. Since the seminal famous paper by Gaspard Monge in 1781 [START_REF] Monge | Mémoire sur la théorie des déblais et des remblais[END_REF], there was a revival of interest in mass transportation in the nineties.

In 1987 [START_REF] Brenier | Décomposition polaire et réarrangement monotone des champs de vecteurs[END_REF][START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF], Brenier proved an existence and uniqueness result of optimal transport maps for the cost c(x, y) = |x -y| 2 in IR n and showed that any such optimal transport map is indeed the gradient of a convex function. Since then, people extended the theory to other costs functions in IR n or to other types of spaces (see [START_REF] Villani | Optimal transport, old and new. Grundlehren des mathematischen Wissenschaften[END_REF]). The aim of this paper is to study existence, uniqueness, and regularity of optimal transport maps for costs functions coming for LQ minimization problems in IR n . Let us consider a linear control system of the form ẋ = Ax + Bu (1.1) where the state x belongs to IR n , the control u belongs to IR m and A, B are n×n and n×m matrices respectively. For every initial state x ∈ IR n and every control u ∈ L 2 [0, 1]; IR m , we denote by x(•; x, u) : [0, 1] → IR n the unique solution to the Cauchy problem ẋ(t) = Ax(t) + Bu(t) for a.e. t ∈ [0, 1],

x(0) = x.

(1.2)

Let us in addition consider a quadratic Lagrangian of the form

L(x, u) = 1 2 x, W x + 1 2 u, U u , (1.3)
where W is a n × n symmetric non-negative matrix and U is a m × m symmetric positive definite matrix. The cost c : IR n × IR n → [0, +∞] associated with (1.1) and (1.3) is given by c(x, y) := inf

1 0 L(x(t; x, u), u(t))dt | u ∈ L 2 [0, 1]; IR m s.t. x(1; x, u) = y , (1.4)
where we set c(x, y) = +∞ if there is no u ∈ L 2 [0, 1]; IR m such that x(1; x, u) = y.

We notice that LQ costs as above include as a particular case the Euclidean cost 1/2|x -y| 2 by taking B = I n , U = I n , A = 0, W = 0 . Costs coming from optimal control have already been studied in [START_REF] Agrachev | Optimal transportation under nonholonomic constraints[END_REF] (see also [START_REF] Figalli | Mass Transportation on sub-Riemannian Manifolds[END_REF]). However, this reference does not give regularity results or properties relating the transport map with the gradient of a convex function. Furthermore, the study of a cost that is finite only for points lying in the same leaf of a foliation, like in the non-controllable case is also original.

The structure of the paper is the following. Our results are stated in Section 2. In Section 3, we introduce some preliminaries in optimal transport theory concerned with Kantorovitch duality. In Sections 4 and 5, we provide the proofs of our results. Then, we present some examples in Section 6. Finally in Section 7, we conclude with several remarks about our results.

Main results

Preliminaries on linear systems

Consider system (1.1) and let V be the smallest linear subspace of IR n that contains the image of the operator B : IR m → IR n and is invariant by A, or, more explicitly,

V = Span IR B, AB, A 2 B, . . . , A n-1 B ⊂ IR n , d = dim V.
(2.1)

The well-known Kalman criterion states that system (1.1) is controllable if and only if d; this is stated below with a description of the situation where d < n. This decomposition can be found for instance in [START_REF] Sontag | Mathematical control theory. Deterministic finite-dimensional systems[END_REF]Lemma 3.3.3].

Proposition 2.1. Let x, y be in IR n . There exists a control u ∈ L 2 [0, 1]; IR m such that x(t; x, u) = y if and only if e A x -y ∈ V (i.e. y lies in the affine subspace e A x + V ).

In particular it exists for all x, y if d; the system -or the pair (A, B)-is then called controllable.

If d < n, after a linear change of coordinates in IR n , the control system (1.1) has the following form:

ẋ = ẋ1 ẋ2 = A 1 A 3 0 A 2 x 1 x 2 + B 1 u 0 , (2.2)
where the state x is partitioned into two blocks x 1 and x 2 of dimension d and n -d respectively and

A 1 , A 2 , A 3 , B 1 are d × d, d × n -d, n -d × n -d and d × m matrices respectively, such that the pair (A 1 , B 1 ) is controllable.

The controllable case

The following result of existence, uniqueness and regularity follows easily from the classical theory of optimal mass transportation. Let µ 0 , µ 1 be two compactly supported probability measures on IR n . Assume that µ 0 is absolutely continuous with respect to the Lebesgue measure L n . Then, there is existence and uniqueness of an optimal transport map T : IR n → IR n . Such a map is characterized by the existence of a convex function ϕ : IR n → IR such that

T (x) = E -1 ∇ϕ(x) for a.e. x ∈ IR n . (2.4)
If in addition, µ 0 , µ 1 are associated with probability densities f 0 , f 1 on supp(µ 0 ) and supp(µ 1 ) respectively, with f 0 and f 1 bounded from below and above, and if supp(µ 0 ) is connected and supp(µ 1 ) convex, then T is continuous.

We postpone several remarks concerning Theorem 2.2 to Section 7.

The non-controllable case

Here d < n in (2.1), hence the evolution of the (n -d)-dimensional x 2 is totally fixed by (2.2) and c(x, y) will obviously be infinite for almost all pairs x, y: these such that y 2 = e A 2 x 2 . A result similar to the above theorem requires first there exists at least one transport map with finite cost, i.e.

inf

S ♯ µ 0 =µ 1 cost c (S) := IR n c(x, S(x)) dµ 0 (x), < +∞. (2.5)
Let π 2 : IR n → IR n-d be the projection on the second block in the coordinates of (2.2): π 2 (x) = x 2 (it is also the projection on the quotient IR n /V ).

Theorem 2.3. Let µ 0 and µ 1 be compactly supported probability measures with continuous

densities f 0 , f 1 (µ 0 = f 0 L n , µ 1 = f 1 L n ).
There exists a transport map with finite cost (2.5) if and only if

(π 2 ) ♯ µ 1 = e A 2 • π 2 ♯ µ 0 . (2.6)
If this is satisfied, then there exists a unique optimal transport map T : IR n → IR n . Moreover, if both densities f 0 , f 1 are bounded from below and above on supp(µ 0 ) and supp(µ 1 ) respectively, and if supp(µ 0 ) and supp(µ 1 ) are both convex, then T is continuous.

Again, remarks concerning Theorem 2.3 are postponed to Section 7.

Preliminaries in optimal transport theory

Given two probability measures µ 0 , µ 1 on IR n and a cost function c :

IR n × IR n → [0, +∞],
we are looking for a transport map T : IR n → IR n which minimizes the transportation cost IR n c(x, T (x)) dµ 0 . The constraint T # µ 0 = µ 1 being highly non-linear, the optimal transport problem is quite difficult from the viewpoint of calculus of variation. The major advance on this problem was due to Kantorovich, who proposed in [START_REF] Kantorovich | On the transfer of masses[END_REF][START_REF] Kantorovich | On a problem of Monge[END_REF] a notion of weak solution of the optimal transport problem. He suggested to look for plans instead of transport maps, that is probability measures γ in IR n × IR n whose marginals are µ 0 and

µ 1 , that is (π 1 ) ♯ γ = µ 0 and (π 2 ) ♯ γ = µ 1 ,
where π 1 : IR n × IR n → IR n are the canonical projections on the first and second variable respectively. Denoting by Π(µ 0 , µ 1 ) the set of plans, the new minimization problem becomes the following:

C µ 0 , µ 1 = min γ∈Π(µ 0 ,µ 1 ) IR n ×IR n c(x, y) dγ(x, y) . (3.1)
If γ is a minimizer for the Kantorovich formulation, we say that it is an optimal plan. Due to the linearity of the constraint γ ∈ Π(µ 0 , µ 1 ), it is simple, using weak topologies, to prove existence of solutions to (3.1) as soon as c is lower semi-continuous (see for instance [START_REF] Villani | Optimal transport, old and new. Grundlehren des mathematischen Wissenschaften[END_REF]). The connection between the formulation of Kantorovich and that of Monge can be seen by noticing that any transport map T induces the plan defined by (Id × T ) ♯ µ 0 which is concentrated on the graph of T . Thus, the problem of showing existence of optimal transport maps can be reduced to prove that an optimal transport plan is concentrated on a graph. Moreover, if one can show that any optimal plan is concentrated on a graph, since γ 1 +γ 2 2 is optimal if so are γ 1 and γ 2 , uniqueness of the transport map easily follows. The following definition is exactly [15, Definition 5.2]:

Definition 3.1. Let c : IR n × IR n → [0, ∞]. A function ψ : IR n → IR ∪ {+∞} is said to be c-convex if it is not identically +∞ and there exists a function ζ : IR n → IR ∪ {-∞, +∞} such that ψ(x) = sup y∈IR n ζ(y) -c(x, y) ∀x ∈ IR n .
Then its c-transform is the function ψ c defined by

ψ c (y) = inf x∈IR n ψ(x) + c(x, y) ∀y ∈ IR n ,
and its c-subdifferential is the set defined by

∂ c ψ := (x, y) ∈ IR n × IR n | ψ c (y) -ψ(x) = c(x, y) .
Moreover, the c-subdifferential of ψ at x ∈ IR n is

∂ c ψ(x) := y ∈ IR n | (x, y) ∈ ∂ c ψ , or equivalently ψ(x) + c(x, y) ≤ ψ(z) + c(z, y) ∀z ∈ IR n .
The functions ψ and ψ c are said to be c-conjugate.

Let us denote by P c (IR n ) the set of compactly supported probability measures. From now on, supp(µ 0 ) and supp(µ 1 ) will denote the supports of µ 0 and µ 1 respectively, i.e. the smallest closed sets on which µ 0 and µ 1 are respectively concentrated. Kantorovitch duality can be stated as follows (see [START_REF] Villani | Optimal transport, old and new. Grundlehren des mathematischen Wissenschaften[END_REF]Theorem 5.10]): Theorem 3.2. Let µ 0 , µ 1 ∈ P c (IR n ) and c : IR n × IR n → [0, ∞) be a lower continuous function. Then there exists a c-convex function ψ : IR n → IR such that the following holds: a transport plan γ ∈ Π(µ 0 , µ 1 ) is optimal if and only if γ(∂ c ψ) = 1 (that is, γ is concentrated on the c-subdifferential of ψ). Moreover ψ can be chosen such that

ψ(x) = sup y∈supp(µ 1 ) {ψ c (y) -c(x, y)} ∀x ∈ IR n , (3.2)
ψ c (y) = inf x∈supp(µ 0 ) {ψ(x) + c(x, y)} ∀y ∈ IR n . (3.3)
By the above theorem we see that, in order to prove existence and uniqueness of optimal transport maps, it suffices to prove that there exist two Borel sets

Z 0 , Z 1 ⊂ IR n , with µ 0 (Z 0 ) = µ 1 (Z 1 ) = 1, such that and ∂ c ψ is a graph inside Z 0 × Z 1 (or equivalently that ∂ c ψ(x) ∩ Z 1 is a singleton for all x ∈ Z 0 ).

Proof of Theorem 2.2

First, we prove that (2.3) holds. This is not new, but the formulas are not given in textbooks in extenso. The case W = 0 is more classical and one may find, for instance in [START_REF] Lee | Foundations of optimal control theory[END_REF] the expression of c in that case:

c(x, y) = x -e -A y , G -1 (x -e -A y) with G = 1 0 e τ A * BU -1 B * e τ A dτ (4.1)
where the controllability Grammian G is positive definite because the system is controllable. This is indeed of the form (2.3). Let us derive the general form from the classical linear Pontryagin Maximum principle (see for instance the same reference).

The pseudo Hamiltonian H 0 : IR n × IR n × IR m → IR associated with the optimization problem under study is given by

H 0 (x, p, u) := p, Ax + Bu - 1 2 x, W x - 1 2 u, U u , (4.2)
the control u that maximizes this expression for each (p, x) is given by

u = U -1 B * p (4.3)
and the Hamiltonian H :

IR n × IR n → IR, defined as H(x, p) = max {H 0 (x, p, u) | u ∈ IR m } is given by H(x, p) = p, Ax - 1 2 x, W x + 1 2 p, BU -1 B * p . (4.4)
Therefore the Hamiltonian differential equation ẋ = ∂H/∂p, ṗ = -∂H/∂x associated to our minimization problem is given by

ẋ ṗ = A BU -1 B * W -A * x p (4.5)
Denote by R(•) : [0, 1] → M 2n (IR) the fundamental solution to the Cauchy problem

Ṙ(t) = A BU -1 B * W -A * R(t) ∀t ∈ [0, 1], R(0) = I 2n , (4.6)
and write it as

R(t) = R 1 (t) R 2 (t) R 3 (t) R 4 (t) ∀t ∈ [0, 1],
where R i (•) is a n × n matrix for i = 1, • • • 4. For every x ∈ IR n fixed, denote by exp x the mapping from IR n to IR n which sends via the Hamiltonian vector field the initial adjoint vector p ∈ IR n to the final state x(1), that is

exp x (p) := R 1 (1)x + R 2 (1)p ∀p ∈ IR n . (4.7)
Controllability of the system (1.1) implies that the affine mapping exp x is onto, and hence a bijection. Indeed, there is, for any x, y ∈ IR n , (at least) one control u ∈ L 2 [0, 1]; IR m which minimizes the cost

1 0 L(x(t; x, u), u(t))dt
among all controls steering x to y in time 1; thanks to the linear maximum principle, there corresponds to each minimizing control a p ∈ IR n such that exp x (p) = y. Hence the matrix R 2 (t) is invertible for all t. Given x and y in IR n , the cost between x and y is given by

c(x, y) = 1 0 1 2 x(t), W x(t) + 1 2 p(t), BU -1 B * p(t) dt,
where x(•), p(•) : [0, 1] → IR n × IR n is the solution to the Hamiltonian system (4.5) satisfying x(0) = x and exp x p(0) = y.

Set p := p(0), that is p = R 2 (1) -1 y -R 1 (1)x . (4.8)
We deduce easily that c is a smooth function of the form

c(x, y) = 1 2 x, Q 1 x + x, Cp + 1 2 p, Q 2 p ∀x, y ∈ IR n , (4.9)
where both Q 1 , Q 2 are non-negative symmetric n × n matrices given by

Q 1 := 1 0 R 1 (t) * W R 1 (t) + R 3 (t) * BU -1 B * R 3 (t) dt, Q 2 := 1 0 R 2 (t) * W R 2 (t) + R 4 (t) * BU -1 B * R 4 (t) dt,
and where C is the n × n matrix given by By definition of the cost c, there holds

C := 1 0 R 1 (t) * W R 2 (t) + R 3 (t) * BU -1 B * R 4 (t) dt.
1 0 L(y(t; y, v), v(t))dt -c y(0; y, v), y ≥ 0 ∀v ∈ L 2 [0, 1], IR m .
Moreover, there is equality in the above inequality whenever v = u. Thanks to the linear maximum principle, since the cost c is smooth, this means that y = exp x (p) with p = -∇ x c(x, y).

Computing ∇ x c(x, y) by differentiating (4.9)-(4.8) with respect to the variable x yields

Q 1 x + Cp -R 1 (1) * R 2 (1) -1 * C * x + Q 2 p = -p,
and finally (recall that Q 1 is symmetric) :

C = -I n + R 1 (1) * R 2 (1) -1 * Q 2 , Q 1 = R 1 (1) * R 2 (1) -1 * C * = CR 2 (1) -1 R 1 (1). (4.10)
Plugging this and (4.8) into (4.9) yields the expression (2.3) for c(x, y) with

D = R 2 (1) -1 R 1 (1) , E = R 2 (1) -1 , F = (R 2 (1) -1 ) * Q 2 R 2 (1) -1 , (4.11)
where D is symmetric because (4.10) implies

Q 1 = -E+E * Q 2 E and Q 1 , Q 2 are symmetric;
Positive definiteness of D and F may be deduced from the fact that there is no solution of (1.1) driving 0 to a nonzero y or a nonzero x to 0 in finite time with a zero cost.

Let us now show the existence and uniqueness of an optimal transport map. Let µ 0 , µ 1 be two compactly supported probability measures in IR n and assume that µ 0 is absolutely continuous with respect to the Lebesgue measure. Let ψ, φ := ψ c : IR n → IR be the Kantorovitch potentials given by Theorem 3.2 (note that c is non-negative valued). First, since c is smooth and both sets supp(µ 0 ), supp(µ 1 ) are assumed to be compact, (3.2)-(3.3) imply that both potentials ψ, φ are locally Lipschitz. Therefore, thanks to Rademacher's Theorem and the fact that µ 0 is absolutely continuous with respect to the Lebesgue measure, ψ is differentiable µ 0 -almost everywhere. Let x ∈ supp(µ 0 ) be a differentiability point of ψ. Let y ∈ supp(µ 1 ) be such that φ(y) -ψ(x) = c(x, y),

and u ∈ L 2 [0, 1], IR m be a control minimizing c(x, y) = inf 1 0 L(x(t; x, v), v(t))dt | v ∈ L 2 [0, 1]; IR m s.t. x(1; x, v) = y .
Then we have

1 0 L(y(t; y, v), v(t))dt ≥ c y(0; y, v), y ≥ φ(y) -ψ y(0; y, v) ∀v ∈ L 2 [0, 1]; IR m , with equality if v = u.
As above, this implies that

y = exp x (p) with p = ∇ψ(x).
This shows that y is uniquely determined for every differentiability point of φ and in turn proves existence and uniqueness of an optimal measurable transport map x → y. Note that the potential ψ can be written as follows:

ψ(x) = sup y∈supp(µ 1 ) {φ(y) -c(x, y)} = - 1 2 x, D x + sup y∈supp(µ 1 )
x, E y -1 2 y, F y + φ(y) .

This shows that the function ϕ : IR n → IR defined by

ϕ(x) := ψ(x) + 1 2 x, D x ∀x ∈ IR n ,
is convex (as the sup of a family of convex functions) while, for almost every x ∈ supp(µ 0 ), deriving the expression of exp x (p) from (4.7) and (4.11),

T (x) = exp x ∇ψ(x) = E -1 Dx + ∇ψ(x) = E -1 ∇ϕ(x) .
It remains to show that under additional assumptions, T is continuous. One obviously has ∇ϕ(x) = E T (x) for all x; it is clear that this and T ♯ µ 0 = µ 1 imply (∇ϕ) ♯ µ 0 = μ1 for the measure μ1 defined by

μ1 (A) = µ 1 E -1 x | x ∈ A
for every Borel set A ⊂ IR n (μ 1 is the pushforward of µ 1 by the map x → Ex). Since ϕ is convex, ∇ϕ is then the optimal transport map from µ 0 to μ1 for the cost ĉ(x, y) = 1 2 |x-y| 2 . This implies the continuity of ∇ϕ -hence of T -according to [START_REF] Villani | Optimal transport, old and new. Grundlehren des mathematischen Wissenschaften[END_REF]Theorem 12.50 (i)] which asserts the following: if μ0 , μ1 are two compactly supported probability measures in IR n associated with densities f0 , f1 which are bounded from below and above on supp( f0 ) and supp( f1 ) respectively, and if supp( f0 ) is connected and supp( f1 ) is convex, then the optimal transport map T : IR n → IR n from μ0 to with respect to the Euclidean quadratic cost ĉ(x, y) = 1 2 |x -y| 2 is continuous. Take μ1 defined above and set μ0 = µ 0 ; they do satisfy these assumptions for µ 0 and µ 1 and applying an invertible linear map does not change convexity of the support of a density, and we just proved that T is ∇ϕ with the same ϕ as above.

Proof of Theorem 2.3

Let us first treat the case d = 0 separately. According to (2.1), either there is no control (m = 0) or the matrix B is zero; in both cases the system reads ẋ = Ax and c(x, y) = Then there is only one map S such that cost c (S) < ∞, given by S(x) = e A x for all x, the compatibility on measures is µ 1 = S ♯ µ 0 with this precise S, and the result is proved in the case d = 0. From now on, we assume that

d ∈ {1, • • • , n -1}, i.e. none of the two blocks in (2.2) is void (d ≥ 1, n-d ≥ 1). For every x = x 1 , x 2 ∈ IR d ×IR n-d = IR n and u ∈ L 2 [0, 1]; IR m , the solution x(•; x, u) : [0, 1] → IR n to the Cauchy problem (1.
2) with the decomposition given by (2.2) satisfies x(t; x, u) = x 1 (t; x, u) x 2 t; x 2 = e tA 1 x 1 + e tA 1 t 0 e -sA 1 A 3 e sA 2 x 2 + B 1 u(s) ds e tA 2 x 2 , (5.1)

for every t ∈ [0, 1]. Denote by x1 (•) = x1 •; x 1 , u : [0, 1] → IR d the solution to the Cauchy problem ẋ1 (t) = A 1 x1 (t) + B 1 u(t) for a.e. t ∈ [0, 1], x1 (0) = x 1 . (5.2)
Then we have

x 1 t; x, u = x1 t; x 1 , u + G(t)x 2 ∀t ∈ [0, 1],
where

G : [0, 1] → M d,n-d (IR) is defined by G(t) := e tA 1 t 0 e -sA 1 A 3 e sA 2 ds ∀t ∈ [0, 1].
Therefore we have for almost every t ∈ [0, 1],

L x(t; x, u), u(t) = 1 2 x(t; x, u), W x(t; x, u) + 1 2 u(t), U u(t) = 1 2 x1 t; x 1 , u , W 1 x1 t; x 1 , u + 1 2 u(t), U u(t) + x1 t; x 1 , u , X t; x 2 + l t; x 2 ,
where the symmetric matrix W is decomposed, in the coordinates of (2.2) as

W = W 1 W 3 W 3 * W 2
with W 1 and W 2 positive definite d × d and n -d × n -d matrices respectively, W 3 a d × n -d matrix, and

X t; x 2 = W 1 G(t)x 2 + W 3 x 2 t; x 2 , l t; x 2 = 1 2
G(t)x 2 x 2 t; x 2 , W 1 W 3 W 3 * W 2 G(t)x 2 x 2 t; x 2 .
In conclusion, we have for every

x = x 1 , x 2 , y = y 1 , y 2 ∈ IR n , c(x, y) = cx 2 x 1 , y 1 -G(1)x 2 + 1 0 l t; x 2 dt if y 2 = e A 2 x 2 +∞ otherwise, (5.3)
where the cost cx 2 :

IR d × IR d → [0, +∞) is defined by cx 2 x 1 , z 1 := inf 1 0 Lx 2 t, x1 (t; x 1 , u), u(t) dt | u ∈ L 2 [0, 1]; IR m s.t. x1 1; x 1 , u = z 1 , (5.4)
for any x 1 , z 1 ∈ IR d , and where the Lagrangian Lx

2 : IR × IR d × IR m → [0, +∞) is defined by Lx 2 t, z, u := 1 2 z, W 1 z + 1 2 u(t), U u(t) + z, X t; x 2 . (5.5)
We proceed now as in the proof of Theorem 2.2. Let x 2 ∈ IR n-d be fixed, the pseudo Hamiltonian H 0 : IR × IR d × IR d × IR m → IR associated with the cost cx 2 is given by

H 0 (t, z, p, u) := p, A 1 z + B 1 u - 1 2 z, W 1 z - 1 2 u, U u -z, X t; x 2 . (5.6)
Then ∂H 0 ∂u = 0 yields u = U -1 B * 1 p and the Hamiltonian H : IR × IR d × IR d → IR is given by

H(t, z, p) = max {H 0 (t, z, p, u) | u ∈ IR m } = p, A 1 z - 1 2 z, W 1 z + 1 2 p, B 1 U -1 B * 1 p -z, X t; x 2 .
Therefore the Hamiltonian system associated to our minimization problem is given by

ż = A 1 z + B 1 U -1 B * 1 p ṗ = -A * 1 p + W 1 z + X t; x 2 .
(5.7)

Denote by R(•) : [0, 1] → M 2d (IR) the fundamental solution to the Cauchy problem

Ṙ(t) = A 1 B 1 U -1 B * 1 W 1 -A * 1 R(t) ∀t ∈ [0, 1], R(0) = I 2d , (5.8)
and write it as

R(t) = R 1 (t) R 2 (t) R 3 (t) R 4 (t) ∀t ∈ [0, 1],
where

R i (•) is a d × d matrix for i = 1, • • • 4. Any solution z(•), p(•) : [0, 1] → IR d × IR d of (5.7
) with z(0) = z, p(0) = p can be written as

z(t) = R 1 (t)z + R 2 (t)p + zx 2 (t) p(t) = R 3 (t)z + R 4 (t)p + px 2 (t) ∀t ∈ [0, 1],
where zx 2 (•), px 2 (•) : [0, 1] → IR d × IR d of (5.7) satisfying zx 2 (0) = 0, px 2 (0) = 0. For every x 1 ∈ IR d , the mapping exp x 1 : IR d → IR d defined by

exp x 1 (p) := R 1 (1)x 1 + R 2 (1)p + zx 2 (1) ∀p ∈ IR d ,
is an affine bijection. For every x 1 , z 1 , p ∈ IR d with z 1 = exp x 1 (p), there holds

cx 2 x 1 , z 1 = 1 2 x 1 , Q 1 x 1 + 1 2 p, Q 2 p + x 1 , C p + x 1 , v x 2 + p, w x 2 + k x 2 , (5.9)
where both Q 1 , Q 2 are non-negative symmetric d × d matrices given by

Q 1 := 1 0 R 1 (t) * W 1 R 1 (t) + R 3 (t) * B 1 U -1 B * 1 R 3 (t) dt, Q 2 := 1 0 R 2 (t) * W 1 R 2 (t) + R 4 (t) * B 1 U -1 B * 1 R 4 (t) dt,
where C is the d × d matrix given by

C := 1 0 R 1 (t) * W 1 R 2 (t) + R 3 (t) * B 1 U -1 B * 1 R 4 (t) dt,
and the vectors v x 2 , w x 2 are given by

v x 2 := 1 0 R 1 (t) * W 1 zx 2 (t) + R 3 (t) * B 1 U -1 B * 1 px 2 (t) + R 1 (t) * X t; x 2 dt, w x 2 := 1 0 R 2 (t) * W 1 zx 2 (t) + R 4 (t) * B 1 U -1 B * 1 px 2 (t) + R 2 (t) * X t; x 2 dt, and
k x 2 := 1 0 1 2 zx 2 (t), W 1 zx 2 (t) + 1 2 px 2 (t), B 1 U -1 B * 1 px 2 (t) + zx 2 (t), X t; x 2 dt.
We proceed now as in the proof of Theorem 2.2. Since cx 2 is smooth, using the linear maximum principle, taking the derivative in the x 1 variable in (5.9) and using that p = R 2 (1) -1 z 1 -R 1 (1)x 1 -zx 2 (1) yields

Q 1 x 1 + C p + v x 2 -R 1 (1) * R 2 (1) -1 * (Q 2 p + C * x 1 + w x 2 ) = -p
and in turn, from (5.3), this implies, for almost every x 2 ∈ IR n-d , that S(x 1 , x 2 ) ∈ N x 2 for almost all x 1 , which means that S(x 1 , x 2 ) has the form (5.12). Now S ♯ µ 0 = µ 1 means R n h(x)dµ 1 (x) = R n h S(x) dµ 0 (x) for any h ∈ L 1 IR n . From Fubini's Theorem together with a change of variable y 2 = e A 2 x 2 in the left-hand side, this yields, for any positive h ∈ L 1 IR n

IR n-d IR d h(x 1 , x 2 )f 1 (x 1 , x 2 )dx 1 dx 2 = IR n-d IR d h S e -A 2 y 2 (x 1 ), y 2 f 0 x 1 , e -A 2 y 2 det e -A 2 dx 1 dy 2 .
This is equivalent to (S x 2 ) ♯ µ x 2 0 = µ x 2 1 and proves the fact. If S satisfies S ♯ µ 0 = µ 1 and cost c (S) < +∞, the form (5.12) implies that, for all measurable h that depends on x 2 only, one has

IR d ×IR n-d h(x 2 )dµ 1 (x 1 , x 2 ) = IR d ×IR n-d h(e A 2 x 2 )dµ 0 (x 1 , x 2 ),
and this implies (2.6). Conversely, assume now that the measures satisfy (2.6); since

x 2 → Mx 2 dµ x 2 0 and x 2 → Nx 2 dµ x 2 1
are the densities of (e A 2 • π 2 ) ♯ µ 0 and (π 2 ) ♯ µ 1 respectively, the measures µ x 2 0 and µ x 2 1 have, for almost all x 2 , the same total mass, hence can be considered as probability measures after normalisation. Then using (5.11) and arguing as in the proof of Theorem 2.2, we deduce that for almost every x 2 ∈ IR n-d , there is a unique optimal transport map T x 2 : M x 2 → N x 2 together with a convex function ϕ x 2 : IR d → IR such that

T x 2 x 1 = E -1 ∇ϕ x 2 (x 1 ), e A 2 x 2 for a.e. x 1 ∈ IR d .
This shows the existence and uniqueness of an optimal transport map T , defined by T (x 1 , x 2 ) = (E -1 ∇ϕ x 2 (x 1 ), e A 2 x 2 ) from µ 0 to µ 1 (and a fortiori shows (2.5)). Note that joint measurability of T with respect to the two variables may be seen as a consequence of the necessity part of [2, Theorem 3.2] which asserts that any optimal plan is concentrated on a c-monotone Borel subset of IR n × IR n . If we assume in addition that f 0 , f 1 are both continuous and bounded from below and above on their support, then each T x 2 is continuous. This is a consequence of Theorem 1 applied for each x 2 to the mass transportation problem from µ x 2 0 to µ x 2 1 with respect to the cost given by (5.11). But by [START_REF] Villani | Optimal transport, old and new. Grundlehren des mathematischen Wissenschaften[END_REF]Corollary 5.23], there is stability of transport maps. If {x k 2 } is a sequence converging to x 2 , the measures µ

x k 2 0 and µ x k 2
1 weakly converge to µ x 2 0 and µ x 2 1 respectively, and if in addition the corresponding costs converge uniformly as well, then the transport {T x k 2 } maps converge to T x 2 in probability. Since all T k x 2 are indeed uniformly Hölder continuous (see [START_REF] Villani | Topics in Mass Transportation[END_REF]Theorem 50 (i)] or [START_REF] Villani | Optimal transport, old and new. Grundlehren des mathematischen Wissenschaften[END_REF]Theorem 12.50]), and references therein), this concludes the proof of Theorem 2.3.

Examples

The case W = 0

If W = 0 then the second line in (4.5) is an independent first order differential equation in the p variable; hence R 3 (t) = 0 for any t ∈ [0, 1] and the form (4.1) of the cost c. If in addition, we assume that A = 0 and that the system is controllable, then the matrix B is necessarily invertible. In that case, we leave the reader to show that the 2n × 2n matrix R(t) has the form

R(t) = I n BU -1 B * 0 n I n .
Then there holds c(x, y) = 1 2 BU -1 B * -1 y -x , y -x ∀x, y ∈ IR n .

And any optimal transport map given by Theorem 2.2 has the form

T = BU -1 B * ∇ϕ,
where we used (4.11), with ϕ is a convex function. This can also be viewed as a consequence of [START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF] because the above cost is the Euclidean norm corresponding to the positive matrix (BU -1 B * ) -1 and is BU -1 B * ∇ϕ is nothing but the gradient of ϕ for this Euclidean metric.

6.2

The case A = 0, B = U = I n

To have a nice form of R(t) let us assume that W is symmetric definite positive. In this case, we leave the reader to show that the 2n × 2n matrix R(t) has the form : R(t) = cosh tW

1 2
sinh tW And any optimal transport map, given by Theorem 2.2, has the form

T = sinh W 1 2 .W -1 2
∇ϕ,

where we used (4.11) and that ϕ is a convex function.

Final comments

It is worth noticing that the existence and uniqueness part in Theorem 2.2 is not new. It is a consequence of [1, Theorem 4.1] together with the Lipschitz regularity of the cost. The formula (2.3) implies that the Ma-Trudinger-Wang tensor associated with the cost is identically zero. Such a result has been obtained previously by McCann and Lee without computing explicitly the cost (see [START_REF] Lee | The Ma-Trudinger-Wang curvature for natural mechanical actions[END_REF]Theorem 1.1]). We finally observe that (2.4) means that the optimal transport map T is convex up to a linear change of coordinates. Such a property is related to the vanishing of the Ma-Trudinger-Wang tensor that we mentioned above and [START_REF] Figalli | Continuity and injectivity of optimal maps for non-negatively cross-curved costs[END_REF]Theorem 4.3]. We refer the interested reader to [START_REF] Figalli | Necessary and sufficient conditions for continuity of optimal transport maps on Riemannian manifolds[END_REF][START_REF] Villani | Optimal transport, old and new. Grundlehren des mathematischen Wissenschaften[END_REF] for a more details about the Ma-Trudinger-Wang tensor and its link with regularity of optimal transport maps.

  e. µ 1 (B) = µ 0 T -1 (B) for all B ⊂ Y measurable), and which minimizes the transportation cost cost c (T ) := X c(x, T (x)) dµ 0 (x).

  Let x, y ∈ IR n and u ∈ L 2 [0, 1]; IR m a control minimizing c(x, y) be fixed. For every control v ∈ L 2 [0, 1], IR m , denote by y(•; y, v) : [0, 1] → IR n the unique solution to the Cauchy problem ẏ(t) = Ay(t) + Bv(t) for a.e. t ∈ [0, 1], y(1) = y.

1 2 1 0

 1 e tA x, W e tA x dt if y = e A x , +∞ otherwise.

and finally Q 1 -R 1 (1) * (R 2 (1) -1 ) * C * = C -R 1 (1) * (R 2 (1) -1 ) * Q 2 +I = v x 2 -R 1 (1) * w x 2 = 0, whence (recall that Q 2 is symmetric):

In conclusion, setting

, and

which in turn implies (by (5.3))

D and F are symmetric definite positive for the same reason as in the proof of Theorem 2.2.

We are now ready to prove the result. For every x2 ∈ IR n-d , set

, and let µ 0 = f 0 L n , µ 1 = f 1 L n with densities f 0 , f 1 ∈ L 1 (IR n ) be two compactly supported probability measures.

Fact: A measurable map S satisfies S ♯ µ 0 = µ 1 and cost c (S) < +∞ if and only if there exists, for almost all x 2 , a measurable S x 2 :

), e A 2 x 2 for a.e. x 1 (5.12)

and S x 2 pushes forward the measure µ x 2 0 on M x 2 defined by

Indeed, by Fubini's Theorem, the map x 1 → S(x 1 , x 2 ) is measurable for almost all x 2 and cost c (S) < +∞ implies IR d c (x 1 , x 2 ), S(x 1 , x 2 ) f 0 (x 1 , x 2 )dx 1 < +∞ for a.e. x 2


Marion Gilson
email: [marion.gilson@cran.uhp-nancy.fr

Hugues Garnier
email: hugues.garnier]@cran.uhp-nancy.fr

Peter C Young
email: p.young@lancaster.ac.uk

Paul M J Van Den Hof

Optimal instrumental variable method for closed-loop identification *

Keywords: System identification, closed-loop identification, optimal instrumental variable, Box-Jenkins models

This paper presents in a new unified way, optimal instrumental variable methods for identifying discrete-time transfer function models when the system operates in closed-loop.

The conditions for the optimal design of prefilters and instruments depending on common model structures are analyzed and different approaches are developed according to whether the controller is known or not. The performance of the proposed approaches is evaluated by Monte-Carlo analysis in comparison with other alternative closed-loop estimation methods.

Introduction

Feedback is present in a variety of practical situations due to safety and/or economic restrictions.

In the last two decades, various attempts have been made to handle linear system identification in the presence of feedback. The main difficulties are well-known and due to the correlation between the disturbances and the control signal induced by the loop. Identification methods dealing with closed-loop experimental data were first developed in the eighties ( [START_REF] Söderström | System identification[END_REF], [START_REF] Ljung | System Identification : Theory for the User -Second Edition[END_REF]); several overviews gather these results (see e.g. [START_REF] Gevers | Towards a joint design of identification and control ?[END_REF], [START_REF] Van Den Hof | Closed-loop issues in system identification[END_REF], [START_REF] Forssell | Closed-loop identification revisited[END_REF]). Renewed attention has been given lately to the problem of closed-loop identification initiated by an emerging interest in identification for model-based control ( [START_REF] Gevers | Identification for control: from the early achievements to the revival of experiment design[END_REF], [START_REF] Landau | Direct controller order reduction by identification in closed loop[END_REF], [START_REF] Goodwin | Bias issues in closed loop identification with application to adaptive control[END_REF]) and/or experiment design in the closed-loop framework ( [START_REF] Mišković | Closed-loop identification of multivariable systems: with or without excitation of all references?[END_REF], [START_REF] Hjalmarsson | From experiment design to closed-loop control[END_REF]). Several methods have been developed recently, therefore, and may be divided into two main approaches depending on the chosen state-space (see e.g. [START_REF] Katayama | Subspace methods for system identification[END_REF] for a survey) or transfer function model form. In this paper, transfer function model identification via instrumental variable (IV) methods is considered.

The IV technique to the identification of transfer function models has indeed a rich history in the control and system literature, with the earliest algorithms of this type dating back to the 1970s for the open-loop case [START_REF] Young | An instrumental variable method for real-time identification of a noisy process[END_REF] and 1980s for the closed-loop situation [START_REF] Söderström | Instrumental variable methods for closed-loop systems[END_REF]. More recently, a so-called 'tailor-made IV algorithm' was proposed [START_REF] Gilson | On the relation between a bias-eliminated least-squares (BELS) and an IV estimator in closed-loop identification[END_REF] where the closed-loop plant is parameterized using the (open-loop) plant parameters. A connexion with the bias eliminated least-squares method was also discussed ( [START_REF] Zheng | Parametric identification of linear systems operating under feedback control[END_REF], [START_REF] Zheng | Application of BELS based methods in direct identification of linear systems from closed loop data[END_REF]). Then, an optimal (minimal) variance result was developed in the closed-loop extended IV identification case, revealing consequences for the choice of weights, filters and instruments. Two bootstrap techniques were developed for estimating AR and ARARX models [START_REF] Gilson | Instrumental variable methods for closed-loop system identification[END_REF].

Attention has been given very recently to more general models and optimal IV solutions are presented either for linear model identification (see e.g. [START_REF] Pintelon | Box-jenkins identification revisited-part i: Theory[END_REF] or [START_REF] Gilson | Instrumental variable methods for closed-loop continuous-time model identification[END_REF] for Box-Jenkins hybrid continuous-time models) or non-linear models [START_REF] Laurain | Refined instrumental variable methods for identifying hammerstein models operating in closed loop[END_REF]. This paper outlines the main aspects of a statistically optimal refined IV approach for the identification of common use structures as ARX, ARARX, OE and BJ, in closed loop. The BJ model structure presents the advantage of allowing the plant and the noise models to have different dynamics and it can be shown that the estimates of the noise model parameters are asymptotically independent of the process model parameter estimates (see [START_REF] Young | The refined instrumental variable method: unified estimation of discrete and continuous-time transfer function models[END_REF] and the prior references therein). An apparent problem with this type of BJ model is that simple IV estimation cannot be used directly because the model is nonlinear-in-the-parameters. However, more sophisticated IV estimation can be used to overcome this limitation [START_REF] Young | The refined instrumental variable method: unified estimation of discrete and continuous-time transfer function models[END_REF]. Indeed, the identification problem is rewritten in this paper to make use of a linear-in-the-parameters predictor next to an additional noise model identification required for determining the optimal prefilter and instrument.

This paper provides, therefore, a survey on the statistically optimal IV solutions for the estimation of discrete-time transfer function models when the system operates in closed-loop, and whatever the controller knowledge is. In this regard, it can be compared to the well-known Prediction Error Minimization (PEM) approach [START_REF] Ljung | System Identification : Theory for the User -Second Edition[END_REF]. However, as illustrated in this paper, the proposed refined IV algorithm seems, as in the open-loop case [START_REF] Young | The refined instrumental variable method: unified estimation of discrete and continuous-time transfer function models[END_REF], less sensitive than the PEM algorithm to the initialization [START_REF] Ljung | Experiments with identification of continuous-time models[END_REF].

The paper is organized as follows. After the preliminaries, the optimal IV solution is presented according to the model structure used, and the controller knowledge. Finally, the different methods are compared with the help of Monte-Carlo simulation examples.

Problem formulation

Cc(q) G 0 (q) H 0 (q) ❣ + - ❣ + + ❣ + + ✲ ❄ ✲ ✻ ✲ ✲ ✲ ✲ ❄ ❄ r 1 (t) r 2 (t) u(t) ξ 0 (t) y(t)
e 0 (t)

Figure 1: Closed-loop system configuration

Preliminaries

Consider a stable, linear, Single Input Single Output (SISO), closed-loop system of the form shown in Figure 1. The data generating system is assumed to be given by the following relations

S :    y(t) = G 0 (q)u(t) + H 0 (q)e 0 (t) u(t) = r(t) -C c (q)y(t), where r(t) = r 1 (t) + C c (q)r 2 (t). ( 1
)
The plant is denoted by G 0 (q) = B 0 (q -1 )/F 0 (q -1 ) with the numerator and denominator degree equals to n 0 , the controller is denoted by C c (q) and q -1 is the delay operator with q -i x(t) =

x(t-i). u(t) describes the plant input signal, y(t) the plant output signal. A coloured disturbance ξ 0 (t) = H 0 (q)e 0 (t) is assumed to affect the closed-loop, where e 0 (t) is a white noise, with zero mean and variance σ 2 e 0 . The following general model structure is chosen to model the system (1)

M : y(t) = G(q, ρ)u(t) + H(q, η)ε(t, θ), (2)
where the parameter vector is given as θ T = (ρ T η T ). The parameterized plant model takes then the form,

G : G(q, ρ) = B(q -1 , ρ) F (q -1 , ρ) = b 1 q -1 + • • • + b n q -n 1 + f 1 q -1 + • • • + f n q -n , (3)
where n denotes the plant model order and with the pair (F, B) assumed to be coprime. The plant model parameters are stacked columnwise in the parameter vector

ρ = [f 1 • • • f n b 1 • • • b n ] T ∈ R 2n . (4)
The plant model order n is assumed known or identified from the data and the parameterized noise model is assumed to be in the form of the following ARMA process,

H : H(q, η) = C(q -1 , η) D(q -1 , η) = 1 + c 1 q -1 + • • • + c m q -m 1 + d 1 q -1 + • • • + d m q -m , (5)
where the associated noise model parameters are stacked columnwise in the parameter vector,

η = d 1 • • • d m c 1 • • • c m T ∈ R 2m . (6)
Note that this paper deals with IV-based methods which are known to give consistent plant model parameter estimates of G 0 irrespective of the structure of H 0 . Now consider the relationship between the plant input and output signals in [START_REF] Forssell | Closed-loop identification revisited[END_REF],

y(t) = G 0 (q)u(t) + H 0 (q)e 0 (t). (7)
If the plant G 0 is included into the chosen model set G (G 0 ∈ G), y(t) can be written as

y(t) = ϕ T (t)ρ 0 + v 0 (t), (8)
where ρ 0 denotes the true plant parameter vector,

ϕ T (t) = [-y(t -1) • • • -y(t -n) u(t -1) • • • u(t -n)] (9)
and v 0 (t) = F 0 (q -1 )H 0 (q)e 0 (t).

In the following, the closed-loop system is assumed to be asymptotically stable and r(t) is an external signal that is persistently exciting of sufficient high order.

The objective is then to estimate the parameter vector from the collected sampled data y(t), u(t) and r(t) when the controller C c (q) is known or not.

Let us introduce some notations used in the following Sections. The closed-loop system can be split up into two different parts as follows

y(t) = G 0 (q) 1 + C c (q)G 0 (q) r(t) + 1 1 + C c (q)G 0 (q) ξ o (t) = ẙ(t) + ỹ(t) u(t) = 1 1 + C c (q)G 0 (q) r(t) - C c (q) 1 + C c (q)G 0 (q) ξ o (t) = ů(t) + ũ(t), (10)
where ẙ(t) and ů(t) are the noise-free input/output, ỹ(t) and ũ(t) the noise part of y(t) and u(t) respectively. The noise-free regressor φT (t) is then defined on the basis of these signals as

follows φT (t) = [-ẙ(t -1) • • • -ẙ(t -n) ů(t -1) • • • ů(t -n)]. (11)
3 IV method for closed-loop identification

Extended IV

The extended-IV estimate is given by (see e.g. [START_REF] Söderström | Instrumental Variable Methods for System Identification[END_REF])

ρxiv (N ) = arg min ρ 1 N N t=1 L(q)ζ(t)L(q)ϕ T (t) ρ - 1 N N t=1 L(q)ζ(t)L(q)y(t) 2 W , (12)
where

ζ(t) ∈ R n ζ with n ζ ≥ 2n is the instrument vector, x 2 W = x T W x,
with W a positive definite weighting matrix and L(q) a stable prefilter.

By definition, when G 0 ∈ G, the extended-IV estimate is consistent under the following two

conditions 1 • ĒL(q)ζ(t)L(q)ϕ T (t) is full column rank, • ĒL(q)ζ(t)L(q)v 0 (t) = 0.

Lower bound for an IV method

The choice of the instrumental variable vector ζ(t), the number of instruments n ζ , the weighting matrix W and the prefilter L(q) may have a considerable effect on the covariance matrix P xiv produced by the IV estimation algorithm. In the open-loop situation the lower bound of the covariance matrix for any unbiased identification method is given by the Cramer-Rao bound (see e.g. [START_REF] Söderström | Instrumental Variable Methods for System Identification[END_REF] and [START_REF] Ljung | System Identification : Theory for the User -Second Edition[END_REF]). The closed-loop situation has been investigated more recently in [START_REF] Gilson | Instrumental variable methods for closed-loop system identification[END_REF].

The main results to be used in the following are only briefly recalled here. It was shown that a minimum value of the covariance matrix P xiv as a function of the design variables ζ(t), L(q) and W exists under the restriction that ζ(t) is a causal function of the external signal r(t) only.

In that case P xiv ≥ P opt xiv with

P opt xiv = σ 2 e 0 [ Ēφ f (t)φ T f (t)] -1 , (13) φf
(t) = L opt (q)φ(t), (14)
L opt (q) = 1 F 0 (q -1 )H 0 (q)
, and

ζ(t) = φ(t). ( 15
)
where φ(t) is the noise-free part of ϕ(t) (see Section 2). Using equations ( 12) and ( 13)-( 15), the following IV estimate is optimal

ρopt (N ) = N t=1 ζ f (t)ϕ T f (t) -1 N t=1 ζ f (t)y f (t) (16)
and where the regressor ϕ f (t) = L opt (q)ϕ(t), the output y f (t) = L opt (q)y(t) and the instrument

vector ζ f (t) = L opt (q)ζ(t) are filtered by L opt (q) (15). 1 The notation Ē[.] = limN→∞ 1 N P N t=1 E[.
] is adopted from the prediction error framework of [START_REF] Ljung | System Identification : Theory for the User -Second Edition[END_REF].

It has to be noted that in the IV estimator considered in ( 16), the optimal choice of instruments and prefilter is dependent on unknown system properties, i.e. the plant as well as the noise dynamics. Whereas dependency of plant dynamics could be taken care of by an iterative procedure where the instrument and prefilter are constructed on the basis of a previous plant model estimate ρi-1 . Knowledge of the noise dynamics is generally missing in an extended IV estimator like [START_REF] Ljung | System Identification : Theory for the User -Second Edition[END_REF] as it is not particularly estimated.

Therefore the next step to an optimal IV method should be to extend the estimator [START_REF] Ljung | System Identification : Theory for the User -Second Edition[END_REF] with a procedure to estimate an appropriate noise model, to be used as a basis for constructing the optimal prefilter L opt (q) given in [START_REF] Mišković | Closed-loop identification of multivariable systems: with or without excitation of all references?[END_REF].

4 Optimal IV identification in closed-loop

IV solutions for commonly used model structures in closed-loop

Optimal IV identification is based on [START_REF] Ljung | System Identification : Theory for the User -Second Edition[END_REF]. However, it requires the knowledge or the estimation of the noise model to build the filter. In this perspective, a general prediction error identification step is added to estimate a model for H 0 (q). This requires the choice of a particular parametrization for H(q, θ), possibly in relation to G(q, ρ).

Two common types of model structure2 have been used so far in the literature [START_REF] Gilson | Instrumental variable methods for closed-loop system identification[END_REF].

• ARX model structure:

A(q -1 , ρ)y(t) = B(q -1 , ρ)u(t) + ε(t, ρ).
In this case the noise model

H(q, ρ) = 1 A(q -1 , ρ) (17)
is already available from the plant model estimate, and no additional noise model estimate is required. On the basis of the estimated plant model G(q, ρ) the corresponding filter is given by

L(q, ρ) = 1 A(q -1 , ρ)H(q, ρ) = 1 (18)
ρ is estimated during the plant model identification step.

• ARARX model structure:

A(q -1 , ρ)y(t) = B(q -1 , ρ)u(t) + 1 D(q -1 , η) ε(t, θ).
In this case the corresponding noise model

H(q, η) = 1 A(q -1 , ρ)D(q -1 , η) (19)
is not available from the plant model ρ only as previously, but an additional noise model estimate is required. In [START_REF] Gilson | Instrumental variable methods for closed-loop system identification[END_REF], this is achieved by identifying D(q -1 , η) as an autoregressive model D(q -1 , η)w(t) = ε(t, η) with

w(t) = A(q -1 , ρ)y(t) -B(q -1 , ρ)u(t), (20)
by using a first plant model estimate ρ and by applying a least square estimator η = arg min

η 1 N N t=1 ε(t, η) 2 .
The optimal filter that results from this procedure is therefore

L(q, η) = D(q -1 , η). ( 21
)
Note that these two methods rely on special structures of the assumed noise models, and therefore limit the possibility to reach consistent estimates of H 0 (q) in situations where these structures are not in accordance with the properties of the underlying systems. Therefore, the IV estimators are then extended to other more general OE and BJ model structures.

• BJ model structure:

The more general structure is obtained by choosing a Box-Jenkins (BJ) model

y(t) = B(q -1 , ρ) F (q -1 , ρ) u(t) + C(q -1 , η) D(q -1 , η) ε(t, ρ, η), (22)
The IV estimator ( 16) cannot be directly applied to this non-linear-in-the-parameter model structure (as in [START_REF] Gilson | Instrumental variable methods for closed-loop system identification[END_REF]). Therefore, the solution is to deal with an IV estimator for the process model coupled with a noise model estimation. This approach is stimulated by the previous development of the refined IV method for open-loop system identification (see [START_REF] Young | Some observations on instrumental variable methods of time-series analysis[END_REF], [START_REF] Young | Refined instrumental variable methods of recursive time-series analysis -part i. single input, single output systems[END_REF],

[26] and [START_REF] Young | The refined instrumental variable method: unified estimation of discrete and continuous-time transfer function models[END_REF] for a unified form that deals with both discrete and continuous-time models).

A natural way to extend the IV estimator [START_REF] Ljung | System Identification : Theory for the User -Second Edition[END_REF] with an estimation of the noise model η is to write

v(t) = C(q -1 , η) D(q -1 , η) ε(t, ρ, η) (23)
with v(t) = y(t) -B(q -1 , ρ)/F (q -1 , ρ)u(t) being available as a measured/reconstructed signal once IV estimator [START_REF] Ljung | System Identification : Theory for the User -Second Edition[END_REF] has delivered a plant model ρ (see the algorithm Section 4.3). Estimation of η in the above equation is then undertaken by an ARMA estimation algorithm on the basis of v(t).

The optimal prefilter that results from this identification procedure is then given by

L(q, ρ, η) = D(q -1 , η) F (q -1 , ρ)C(q -1 , η) . ( 24
)
• OE model structure:

When choosing an Output Error (OE) model structure for the noise dynamics, we arrive at

y(t) = B(q -1 , ρ) F (q -1 , ρ) u(t) + ε(t, ρ), (25)
and of course no additional estimation of the noise dynamics is required. The optimal prefilter however changes, and now becomes given by

L(q, ρ) = 1 F (q -1 , ρ) . ( 26
)

Construction of the instruments

According to [START_REF] Mišković | Closed-loop identification of multivariable systems: with or without excitation of all references?[END_REF], the instruments make use of the noise-free regressor φ(t) which is unknown.

It has to be estimated, therefore, from data during the plant model estimation. Two situations have to be investigated according to the controller knowledge.

Controller unknown

When the controller is unknown, the IV series ζ f (t, θ) are generated on the basis of the two closedloop auxiliary models denoted as G yr (q, ρyr ) and G ur (q, ρur ) (estimated TF between y(t)/r(t)

and u(t)/r(t) respectively) as follows ŷ(t) = G yr (q, ρyr )r(t) [START_REF] Zheng | Application of BELS based methods in direct identification of linear systems from closed loop data[END_REF] û(t) = G ur (q, ρur )r(t) (28)

ζ T f (t, θ) = L(q, θ) -ŷ(t -1) . . . -ŷ(t -n) û(t -1) . . . û(t -n) ( 29
)
where the closed-loop parameters ρyr and ρur can be estimated by the following IV technique

ρyr = N t=1 ζ r (t)ϕ T yr (t) -1 N t=1 ζ r (t)y(t) (30) ρur = N t=1 ζ r (t)ϕ T ur (t) -1 N t=1 ζ r (t)u(t) (31)
with the instruments ζ r (t) generated from the excitation signal r(t).

Controller known

When the controller is known, it is worthwhile to use this information into the identification procedure. Therefore, it may be used with the open-loop plant model to construct the instruments and the filter while satisfying the optimal conditions (13)- [START_REF] Mišković | Closed-loop identification of multivariable systems: with or without excitation of all references?[END_REF]. Indeed, the two closed-loop transfer functions (G yr (q, ρyr ) and G ur (q, ρur )) are no longer needed and the noise-free estimation of the instrumental vector can be achieved by using the previous-step estimated open-loop transfer function. In this case, the IV series ζ f (t, θ) are generated on the following basis

ŷ(t, ρ) = G(q, ρ) 1 + C c (q)G(q, ρ) r(t), (32)
û(t, ρ) = 1 1 + C c (q)G(q, ρ) r(t), (33)
ζ T f (t, θ) = L(q, θ) -ŷ(t -1, ρ) . . . -ŷ(t -n, ρ) û(t -1, ρ) . . . û(t -n, ρ) (34)

Algorithm

The outline of the optimal IV algorithm for any model structures described Section 4.1, within the closed-loop context is given below.

Step 1. Initialization

• Estimate a first plant model, get ρ0 and denote G(q, ρ0 ) = B(q -1 , ρ0 )/F (q -1 , ρ0 ) the corresponding transfer function.

• Set the initial noise model estimates C(q -1 , η0 ) = D(q -1 , η0 ) = 1 and i = 1.

• If the controller is unknown, estimate as well the two closed-loop models (30)-(31) used in the following to build the instruments.

Step 2. Estimate by IV Generate the filtered instruments 3 according to the model structure used and the previous estimated models as    L(q, θi-1 ) computed using either ( 18), ( 21), [START_REF] Young | The refined instrumental variable method: unified estimation of discrete and continuous-time transfer function models[END_REF], or (26)

ζ f (t, θi-1 ) = L(q, θi-1 ) -ŷ(t -1) . . . -ŷ(t -n) û(t -1) . . . û(t -n) (35)
where θi-1,T = ρi-1,T ηi-1,T , and ζ f (t, θi-1 ) can be seen as a filtered estimate of the noise-free part of the regressor vector ϕ(t) (9) based on estimates ŷ(t) and û(t) of the noise-free output and input of the plant, respectively (using ( 27)-( 28) if the controller is unknown and (32)-(33

)
if it is known). Determine the IV estimate using the prefilter and these instruments ρi =

N t=1 ζ f (t, θi-1 )ϕ T f (t, θi-1 ) -1 N t=1 ζ f (t, θi-1 )y f (t, θi-1 ) (36)    ϕ f (t, θi-1 ) = L(q, θi-1 )ϕ(t)
y f (t, θi-1 ) = L(q, θi-1 )y(t) (37)

3 i stands for the ith iteration

This yields B(q -1 , ρi ) and F (q -1 , ρi ). Denote the corresponding transfer function by G(q, ρi ) = B(q -1 , ρi )/F (q -1 , ρi ).

Step 3. Obtain an estimate of the noise model parameter vector η i based on the estimated noise sequence

Use one of the noise model identification scheme 4 described in Section 4 to estimate ηi and the associated transfer function H(q, ηi ).

Step 4. Repeat from step 2. Stop when F (q -1 , ρ), B(q -1 , ρ), H(q, η) and L(q -1 , θ) have converged.

Step 5. Compute the estimated parametric error covariance matrix Pθ associated with the parameter estimates, from

Pθ = σ2 ϕ T f (t, θ)ϕ f (t, θ) -1 , ( 38
)
where σ2 is the sample variance of the estimated residuals.

Then, according to the model structure used, the resulting algorithm will be referred to as cliv arx , cliv ararx , cliv bj or cliv oe .

Remark -Initialization step

As illustrated in the following simulation example Section, the proposed optimal IV estimation method is robust to the initialization. Therefore, several methods may be used in this initialization step of the algorithm. Indeed, this first plant estimate is only used to construct the first filter and instruments, and therefore this suggests that minor deviations from the optimal value only causes second-order effects in the resulting accuracy of the first iteration model.

Simulation examples

These examples are concerned with a simulation model based on the following relations [START_REF] Forssell | Closed-loop identification revisited[END_REF], where

G 0 (q) = 0.0997q -1 -0.0902q -2 1 -1.8858q -1 + 0.9048q -2 , n = 2 C c (q) = 10.75 -9.25q -1 1 -q -1 ,
4 One of the advantages of the proposed algorithms is that they provide consistent plant estimates while still exploiting the pseudo-linear regression type of estimation. Indeed, the IV based pseudo-linear regression method recently suggested in [START_REF] Young | An instrumental variable approach to ARMA model identification and estimation[END_REF] could be used to estimate the ARMA process noise in this step of the algorithms.

The excitation signal r(t) is a pseudo random binary signal of maximal length, with the number of stages for the shift register set to 9 and the clock period set to 8; e 0 (t) is a white noise uncorrelated with r(t).

All of the IV methods presented in this paper lead to consistent results if the assumption G 0 ∈ G is fulfilled and optimal estimates if S ∈ M. This example section illustrates this point with the use of the two following examples.

Example 1: white measurement noise

Firstly, a white noise disturbance (H 0 (q) = 1) is considered in order to evaluate the performance of the cliv oe algorithm in the case S ∈ M. The plant parameters are estimated on the basis of closed-loop data of length N = 4088. A Monte-Carlo simulation of 100 runs is used for a signal-to-noise (SNR) ratio given as SNR = 10 log

P ẙ P e = 35 dB, (39)
where P e represents the average power of the zero-mean additive noise on the system output (e.g. the variance) while P ẙ denotes the average power of the noise-free output fluctuations. The proposed method is compared to:

• clivr : the first basic IV method developed to handle the closed-loop case which uses the delayed version of the reference signal as instruments [START_REF] Söderström | Instrumental variable methods for closed-loop systems[END_REF];

• cliv arx : the optimal IV method for an ARX model structure making use of the controller knowledge (see Section 4);

• cliv oe1 : the proposed IV method for OE model, assuming the controller unknown (see Section 4.3);

• cliv oe2 : the proposed IV method for OE model, making use of the controller knowledge (see Section 4.2.2);

• pem: direct use of the closed-loop data (u(t) and y(t)) in the pem (OE here) method (see e.g. [START_REF] Forssell | Closed-loop identification revisited[END_REF]). The method is known to be theoretically efficient in the S ∈ M case.

Monte-Carlo simulation (MCS) results are presented in Table 1 where the mean and standard deviation of the estimated parameters are displayed. The average number of iterations (N iter )

for the pem and cliv oe algorithms are also given. All the methods provide unbiased results with smaller standard deviations for the two proposed cliv oe and pem methods. Indeed, these methods lead to better results thanks to the iterative estimation procedure, even though the number of iterations required for convergence is quite low (although it should be noted that, on average, pem requires more iterations to converge here). As expected, the clivr method provides the least accurate results since it is a basic (and not optimal) IV approach.

Moreover, it has to be noted that, contrary to the cliv oe and pem algorithms, the cliv arx method only works within the G 0 ∈ G assumption (and not S ∈ M). However, even if the standard deviations are higher than those of the cliv oe and pem, this refined IV provides satisfactory results and is robust to a miss-specified noise model.

H 0 (q) = 1 + 0.5q -1 1 -0.85q -1 .
The following algorithms are used to estimate this model:

• the clivr method (with G 0 ∈ G, S ∈ M);

• the cliv ararx method (with G 0 ∈ G, S ∈ M);

• the cliv oe1 method (with G 0 ∈ G, S ∈ M);

• cliv bj1 : the proposed IV method for BJ model, assuming the controller unknown (with S ∈ M);

• cliv bj2 : the proposed IV method for BJ model, making use of the controller knowledge (S ∈ M);

• the pem algorithm (with S ∈ M).

The plant parameters are estimated on the basis of closed-loop data of length N = 4088. A Monte-Carlo simulation of 100 runs is performed for SNR = 25dB. As previously, the results obtained using different algorithms are presented in Table 2 along with the average number of iterations needed for the cliv bj and pem algorithms. The Bode diagrams of the 100 models identified by the clivr, pem, cliv ararx and cliv bj algorithms are displayed in Figure 2.

As previously, all the IV-based methods give unbiased results. However, the basic clivr does clearly not deliver statistically efficient results. The additional noise model estimation in the refined IV type of algorithms allows to achieve statistically accurate results, especially with the cliv bj methods which makes it possible to exactly model the noise. Moreover, it could be noted that the cliv arx is again robust to a miss-specified noise model. Furthermore, in Table 2 and Figure 2, it could be noticed that the pem algorithm (in the situation S ∈ M) is not able to converge to the global minimum at each run and therefore leads sometimes to erroneous results (see the large standard deviation of the pem estimates in Table 2). This is an expected result since it is well-known that this method may suffer from initialization issues (see e.g. [START_REF] Ljung | Initialisation aspects for subspace and output-error identification methods[END_REF] and [START_REF] Ljung | Experiments with identification of continuous-time models[END_REF]). However, better results may be obtained using the pem routine along with a more sophisticated multiple initialization step.

Additionally, it can be seen from Figure 2 that the proposed IV method named cliv oe still gives really accurate results (unbiased and low variance) even in the G 0 ∈ G but S ∈ M case.

Obviously, the variance is a little bit higher than with cliv bj but still satisfactory. This point illustrates the 'robustness' of the method face to a miss-specified noise model.

G 0 ∈ G case

This section illustrates the performance of the proposed method face to miss-specified process and noise model structures (with S ∈ M and G 0 ∈ G). The example is the same as in Section 5.2.1. A single experience is conducted using the cliv oe method, estimating a third order model and assuming a white noise model.

The results are given in Figure 3 where the comparison between the model output for the input signal and the measured output are depicted (over the first 1000 data points). Since we generated the data, we enjoy the luxury of comparing the model output to the noise-free system output (using the compare Matlab routine). The two plots coincide quite well and the coefficient of determination is equal to 98.9% (the pem routine on the same data set and assuming the same model structure leads to 81.4%).

Conclusion

This paper has highlighted the use of optimal instrumental variable techniques in the closed-loop identification framework. Several available techniques have been presented in a new unified way according to the chosen model structure. Furthermore, a method able to estimate efficiently closed-loop Box-Jenkins models has been proposed. The general Box-Jenkins model structure identification has been handled by making use of a linear-in-the-parameters predictor next to an additional noise model identification required for determining the optimal filter and the instruments. The proposed method presents several advantages: it does not need to know the controller, it may be applied whatever the structure (linear or not) of the controller, it is able to handle the estimation of the more common use model structures (as ARX, ARARX, OE, BJ), it is quite robust to the initialization step and it does not suffer from convergence problem.

Figure 2 :

 2 Figure 2: Estimated Bode diagrams (gain and phase (degree)) of the plant model G(q, η) over the 100 MCS, colored noise

Figure 3 :

 3 Figure 3: Estimated (with cliv oe ) and ∈ G case

Table 1 :

 1 Nevertheless, the proposed cliv oe algorithms allow to achieved more accurate results thanks to the appropriate noise model estimation. The two cliv oe methods (with or without the knowledge of the controller) perform similarly overall. Mean and standard deviation of the 100 estimated models, white measurement noise

	parameters	b1	b2	â1	â2	Niter
	true values	0.0997	-0.0902	-1.8838	0.9048
	clivr	0.0997 ± 0.30e-3 -0.0898 ± 3.80e-3	-1.8821 ± 35e-3	0.9014 ± 32.4e-3
	clivarx	0.0997 ± 0.30e-3 -0.0900 ± 1.10e-3	-1.8772 ± 9e-3	0.8977 ± 8.5e-3
	clivoe1	0.0997 ± 0.05e-3 -0.0902 ± 0.06e-3 -1.8858 ± 0.12e-3 0.9048 ± 0.10e-3	3.53
	clivoe2	0.0997 ± 0.05e-3 -0.0902 ± 0.06e-3 -1.8858 ± 0.10e-3 0.9048 ± 0.09e-3	3.89
	pem	0.0996 ± 0.05e-3 -0.0901 ± 0.06e-3 -1.8858 ± 0.12e-3 0.9048 ± 0.10e-3	4.14
	5.2 Example 2: colored measurement noise
	5.2.1 G 0 ∈ G case
	A second example is used to illustrate the performance of the proposed methods in the case of
	a colored measurement noise, with

Table 2 :

 2 Mean and standard deviation of the 100 estimated models, colored measurement noise

According to standard PE identification literature the autoregressive part of the following models is denoted by polynomial A (in contrast with the notation F in section 2).


Séverine Rosat
email: severine.rosat@eost.u-strasbg.fr

Nicolas Florsch

Jacques Hinderer

Muriel Llubes

Estimation of the Free Core Nutation parameters from SG data: sensitivity study and comparative analysis using linearized Least-Squares and Bayesian methods

Keywords: Free Core Nutation, Superconducting Gravimeter, Bayesian analysis

The Free Core Nutation (FCN) is investigated with the help of its resonance effect on the tidal amplitudes in Superconducting Gravimeter (SG) records of the GGP network. The FCN resonance parameters are combined in a resonance equation involving the Earth's interior parameters. The sensitivity of the FCN parameters to the diurnal tidal waves demonstrates that the quality factor of the FCN is strongly dependent on the accuracy of the imaginary part estimates of the gravimetric factors close to the resonance. The weak amplitude of Ψ 1 tidal wave on the Earth, which is the closest in frequency to the FCN, in addition to errors in ocean loading correction, explains the poor determination of the quality factor Q from surface gravimetric data. The inversion of tidal gravimetric factors leads to estimates of the period, Q and resonance strength of the FCN. We show that, by inverting log(Q) instead of Q, the results using the least-squares method optimized using the Levenberg-Marquardt algorithm are in agreement with the Bayesian probabilistic results and agree with the results obtained A c c e p t e d M a n u s c r i p t 2 from VLBI nutation data. Finally, a combined inversion of 7 GGP European SG data is performed giving T = 428 +/-3 days and 7762 < Q < 31989 (90% C.I.

). An experimental estimate of the internal pressure Love number is also proposed.

Introduction

Because of the fluidity of the core, the Earth has a rotational mode, called the free core nutation (FCN) with a period almost diurnal in Earth-fixed coordinates. The FCN parameters (period, damping) strongly depend on the coupling mechanism at the core-mantle boundary (flattening, topography, electro-magnetic coupling…). The FCN can be detected by its effect on the Earth's rotation, using the VLBI network analyses, or by studying its effects on the gravity field. As the tidal potential contains some diurnal components, a resonance occurs in the diurnal frequency band. This resonance effect can be observed in time-varying gravity data continuously recorded on the Earth's surface by Superconducting Gravimeters (SGs) of the Global Geodynamics Project (GGP) network [START_REF] Crossley | Network of superconducting gravimeters benefits a number of disciplines[END_REF]. The FCN resonance in gravity data is commonly represented by a damped harmonic oscillator model that we invert in order to determine the FCN frequency, quality factor Q and the transfer function of the mantle (or the resonance strength). The usual approach to solve this non-linear inverse problem is to use a linearized least-squares method optimized based on the Levenberg-Marquardt algorithm (Marquardt, 1963 -Numerical Recipes Fortran Chapter 15.5-see for instance [START_REF] Defraigne | Stacking gravity tide measurements and nutation observations in order to determine the complex eigenfrequency of the nearly diurnal free wobble[END_REF][START_REF] Defraigne | Correction to 'Stacking gravity tide measurements and nutation observations in order to determine the complex eigenfrequency of the nearly diurnal free wobble[END_REF][START_REF] Sato | Parameters of the fluid core resonance inferred from superconducting gravimeter data[END_REF]Ducarme et al., 2007). However [START_REF] Florsch | Bayesian estimation of the free core nutation parameters from the analysis of precise tidal gravity data[END_REF] have demonstrated the inadequacy of using such a least-squares method, because the statistical distribution of Q is definitely not Gaussian. They have A c c e p t e d M a n u s c r i p t 3 proposed instead the use of a Bayesian approach to invert the FCN parameters, since the Bayesian method better propagates the information to the parameters. [START_REF] Neuberg | Stacking gravity tide observations in Central Europe for the retrieval of the complex eigenfrequency of the nearly diurnal free wobble[END_REF] first proposed an inversion of stacked gravity tide measurements in central Europe to retrieve the FCN parameters using the Marquardt optimized linearized least-squares. Then [START_REF] Defraigne | Stacking gravity tide measurements and nutation observations in order to determine the complex eigenfrequency of the nearly diurnal free wobble[END_REF] extended the gravity stack to the nutation observations. In those past studies, the obtained Q-value was abnormally small and sometimes even negative. [START_REF] Sato | Parameters of the fluid core resonance inferred from superconducting gravimeter data[END_REF] used 1/Q instead of Q as a parameter to be inverted using a modified Marquardt least-squares method since 1/Q seems to be Gaussian [START_REF] Sato | Resonance parameters of nearly diurnal free core nutation measured with three superconducting gravimeters in Japan[END_REF].

However they obtained a Q-value still smaller than the one retrieved from the VLBI nutation analysis (Table 1). The first application of the Bayesian method was proposed by [START_REF] Florsch | Bayesian estimation of the free core nutation parameters from the analysis of precise tidal gravity data[END_REF], who introduced log(Q) as a parameter instead of Q, in order to preserve the positivity of Q, and obtained a Q-value greater than 20,000. More recently, [START_REF] Ducarme | European tidal gravity observations: Comparison with Earth Tides models and estimation of the Free Core Nutation (FCN) parameters[END_REF] inverted log(Q) using both a Bayesian and a least-squares approach but applied on averaged gravimetric factors from European sites. They obtained a value for Q consistent with the VLBI result. The other studies based on the least-squares method are summarized in Table 1. Note that [START_REF] Koot | Estimation of Earth interior parameters from a Bayesian inversion of very long baseline interferometry nutation time series[END_REF] performed an estimation of the FCN resonance parameters from VLBI nutation series using a Bayesian statistical approach in the time domain.

Here we propose a comparison of the results given by the linearized least squares method optimized by the Levenberg-Marquardt algorithm with the Bayesian inversion applied on SG gravity records. We show that the Q-value obtained from SG data is now in agreement with the value inverted from VLBI nutation series whatever the method used, least-squares or Bayesian inversion. Besides, we demonstrate that the poor constraint on the Q-value obtained by [START_REF] Florsch | Bayesian estimation of the free core nutation parameters from the analysis of precise tidal gravity data[END_REF] was due to the large uncertainty on the phase of the diurnal tidal waves close to the resonance. In the first part we describe the FCN resonance model. Then, we review the theory of the Bayesian method and the Levenberg-Marquardt optimization algorithm applied to linearized least-squares. A qualitative study is then performed to check the sensitivity of the gravity factors to the FCN parameters. Finally, we invert the FCN resonance parameters using a combination of 7 European SG time-series.

The FCN resonance model

The basic equation used to describe the resonance of the FCN in the tidal gravity is usually written as (Hinderer et al., 1991a):

nd j ref j a σ σ δ δ - + = , (1)
where j δ ~corresponds to the complex gravimetric factor observed for every tidal wave of frequency σ j , expressing the damping due to all the physical processes involved in the resonance, is defined as

I nd R nd Q σ σ 2 =
. The quantities R nd σ and I nd σ are positive by definition, therefore they should follow a log-normal distribution law [START_REF] Tarantola | Inverse Problem Theory and Model Parameter Estimation[END_REF][START_REF] Florsch | Bayesian estimation of the free core nutation parameters from the analysis of precise tidal gravity data[END_REF] to avoid possible negative values. It is therefore recommended to include the a priori positivity of Q in the model by changing the variable Q = 10 x and inverting for x, instead of Q. [START_REF] Florsch | Bayesian estimation of the free core nutation parameters from the analysis of precise tidal gravity data[END_REF] performed the inversion by treating ref δ ~as an unknown and showed that a correlation exists between the real parts R ref δ and R nd σ , and between a R and R ref δ , but the correlation between a R and T is much stronger. As ref δ ~ has a weak influence on the values of T and Q, we do not include this parameter in the inversion process. In previous studies (e.g. [START_REF] Defraigne | Stacking gravity tide measurements and nutation observations in order to determine the complex eigenfrequency of the nearly diurnal free wobble[END_REF][START_REF] Defraigne | Correction to 'Stacking gravity tide measurements and nutation observations in order to determine the complex eigenfrequency of the nearly diurnal free wobble[END_REF]Ducarme et al., 2007), the observed value for the tidal wave O 1 was used as the reference gravimetric factor. In our case we will use the mean value of the theoretical inelastic amplitude factors of the O 1 and OO 1 waves computed for the Wahr-Dehant model [START_REF] Wahr | Body tides on anelliptical, rotating, elastic and oceanless Earth[END_REF][START_REF] Dehant | Tidal parameters for an inelastic Earth[END_REF]. By doing so, we suppose that the scale factors of the SG used here are accurate enough, which is usually true (better than 0.3% accuracy, e.g. [START_REF] Amalvict | Absolute Gravity Measurements and Calibration of SG-CT031 at Canberra, 1999-2000[END_REF][START_REF] Sato | Parameters of the fluid core resonance inferred from superconducting gravimeter data[END_REF]. We could also have normalized by the observed O1 amplitude as done by [START_REF] Sato | Parameters of the fluid core resonance inferred from superconducting gravimeter data[END_REF] but we suppose that the scaling error is negligible with respect to the ocean loading uncertainty.

The equations relative to the resonance model are written:

A c c e p t e d M a n u s c r i p t

6                      + - - + =         + - - - + = - - - - 2 2 2 2 2 10 ) ( ) ( 2 10 2 10 ) ( 2 10 ) ( x R nd R nd j R nd j I x R nd R I j x R nd R nd j x R nd I R nd j R ref R j a a a a σ σ σ σ σ σ δ σ σ σ σ σ σ δ δ with 2 1 1 R OO R O ref δ δ δ + = .
(2)

Equations ( 2) will be our formulation for the resonance model. [START_REF] Hinderer | Resonance effects of the Earth's fluid core[END_REF], [START_REF] Hinderer | Elasto-gravitational deformation, relative changes in gravity and earth dynamics[END_REF], [START_REF] Defraigne | Détermination de la fréquence et du facteur de qualité du "Nearly-Diurnal Free Wobble" en utilisant conjointement des observations de marées et de nutations[END_REF] and [START_REF] Legros | The influence of the solid inner core on gravity changes and spatial nutations induced by luni-solar tides and surface loading[END_REF] have written the analytical expression for the resonance strength as a function of the Earth's interior parameters:

      - Ω -         -         - - = e e A A A A e A A e A A a m nd m f m γ δ σ γ σ δ 1 ' 1 ~1 2 (3)
where γ is a compliance [START_REF] Mathews | Modeling of nutation and precession: New nutation series for nonrigid Earth and insights into the Earth's interior[END_REF] related to the geodetic parameter q 0 (ratio of centrifugal acceleration to gravitational acceleration) and the Love number h f by

f h q 2 0 = γ . 1
δ is the pressure Love number that represents the gravity response to the pressure potential acting at the CMB (Hinderer et al. 1991a). e is the dynamic ellipticity of the Earth, Ω its rotation rate and A/A m is the ratio of Earth to mantle moments of inertia.

( )

β σ - Ω - = f m nd e A A ' ~
is the FCN frequency in the inertial reference frame ( 1993) also considered the effect of the inner core including the Free Inner Core Nutation (FICN) resonance effect. We neglect it here as the effect would be too small to be seen in SG records.

Ω + = nd nd σ σ ' ~),
In the following, we will estimate the period T, the quality factor Q, as well as the real and imaginary parts of the resonance strength a ~using two methods: a probabilistic approach based on the Bayesian inversion and the more classical linearized least-squares optimized using the Levenberg-Marquardt algorithm.

A review of the methods

a. The Bayesian approach

The Bayesian inversion consists in propagating the information (or knowledge) provided by the measurements through an assumed physical model (perfectly or probabilistically known) to the parameters and to include the a priori knowledge of the model parameters. Both the data and the model parameters are described with probability distributions. The Bayesian approach preserves the full knowledge provided by the data combined with the physical law and the a priori information on the data and model parameters. Therefore, it is the most suitable method to perform the inversion of non-linear problems (Tarantola and Valette, 1982a, b). For more details about the Bayesian method, please refer to [START_REF] Florsch | Bayesian estimation of the free core nutation parameters from the analysis of precise tidal gravity data[END_REF].

The Bayesian probability distribution of the parameter vector θ is given by:

                          ∆ - +         ∆ - - = = ∑ j j j th j j j th j I R R nd k a a x p p 2 2 ) Im( ) Im( ) Im( ) Re( ) Re( ) Re( 2 1 exp ) , , , ( ) ( δ δ δ δ δ δ σ θ , ( 4
)
where k is a normalisation factor in order that the integral of this equation is unity , Re denotes the real part and Im the imaginary part, δ j holds for the j th measurement value of the The previous formula gives the general probability laws for the parameter vector θ. In order to obtain the law for one or two parameters, we compute the marginal pdfs by integration of the probability function over selected parameters. For instance, the joint pdf integrated with respect to R nd σ is defined by:

∫ = R nd I R nd I R a a x d a a x p a a x p I R σ σ ) , , , ( ) , , ( R , , .
Notice that when integrating over selected parameters, we limit a priori the space parameter on a finite domain. Two further integrations of the pdf lead to the marginal probability law for each of the parameters.

b. The linearized least-squares approach

As for the Bayesian method, the linearized least-squares approach is a subset of the non-linear least squares generalized by Tarantola and Valette (1982b). The only difference is how to treat a priori information, and what pdfs are involved. The least-squares method is based on the determination of the best-fit parameters by minimizing the merit function:

∑ ∆         - + - = j j nd j ref j a 2 2 2 ) ( δ σ σ δ δ χ (5)
However, when the model is non-linear, the minimization must proceed iteratively given initial values for the parameters. The procedure finishes when χ 2 stops decreasing. The χ 2 function is linearized to a quadratic form depending on the Hessian matrix for the inverse-Hessian method or the gradient of χ 2 for the steepest descent method. The Levenberg-Marquardt (also called Marquardt method; [START_REF] Marquardt | An algorithm for least-squares estimation of non-linear parameters[END_REF] algorithm is in fact a smooth variant between these two methods. The Marquardt method works very well in practice and A c c e p t e d M a n u s c r i p t 9 has become a standard of non-linear least-squares routines. It has been widely used in previous FCN retrieval studies like in [START_REF] Neuberg | Stacking gravity tide observations in Central Europe for the retrieval of the complex eigenfrequency of the nearly diurnal free wobble[END_REF], [START_REF] Richter | Chandler effects and the nearly diurnal free wobble as determined from observations with a superconducting gravimeter[END_REF], [START_REF] Zürn | Investigation of the nearly diurnal free wobble resonance in individual tidal records[END_REF], and [START_REF] Defraigne | Stacking gravity tide measurements and nutation observations in order to determine the complex eigenfrequency of the nearly diurnal free wobble[END_REF][START_REF] Defraigne | Correction to 'Stacking gravity tide measurements and nutation observations in order to determine the complex eigenfrequency of the nearly diurnal free wobble[END_REF]. We refer to the Numerical Recipes [START_REF] Press | Numerical Recipes in FORTRAN: the Art of Scientific Computing[END_REF] Chapter 15.5 for the optimization subroutines.

The set of linear equations

∑ = ∆ j k j kj β θ α
is solved for the increments ∆θ j that, added to the current trial parameters, give the next approximation. The matrix α is called the curvature matrix and is equal to one-half times the Hessian matrix:

j k kj a a ∂ ∂ ∂ = 2 2 2 1 χ α and k k a ∂ ∂ - = 2 2 1 χ β
, where a k (or a j ) corresponds to one parameter. In the Marquardt optimization algorithm, the increments are related to the curvature matrix through a "damping factor" λ as

j jj j β λα θ 1 = ∆ .
When λ is very large, the diagonal of the curvature matrix is dominant and the solution tends to the initial parameters (the increments tend toward zero). On the other hand, as λ approaches zero, the increments become very large and we may have a divergence problem. The damping factor λ is also known as the Marquardt factor. The damping factor effectively constrains the range of values that the increments ∆θ j can take.

The model used for the inversion is described by equation (2). In order to impose the positivity of the quality factor Q, we perform the change of variable x = log 10 (Q), as for the Bayesian approach, then we minimize equation ( 5). We also estimate the damping factor λ, which best minimizes χ 2 .

Sensitivity and correlations between the FCN parameters

In order to check the sensitivity of the diurnal tidal wave amplitudes to the FCN parameters, we perform some qualitative studies of the resonance transfer function presented in equation ( 2) and based on the previously estimated values. We have s rad deg/h and a I ≈ 2 10 -10 rad/s ≈ 4 10 -5 deg/h (e.g. [START_REF] Hinderer | Elasto-gravitational deformation, relative changes in gravity and earth dynamics[END_REF][START_REF] Florsch | Bayesian estimation of the free core nutation parameters from the analysis of precise tidal gravity data[END_REF]; see next section).

From equation ( 2) and with these values, for any diurnal waves of angular frequency j σ , we

can see that the quantity Q a R nd I 2 σ is negligible with respect to ) ( R nd j R a σ σ - and ( ) 2 2 2 R nd j R nd Q σ σ σ - <<        
, hence the amplitude gravimetric factor can be approximated by:

         - - + ≈ - + ≈ - 2 ) ( ) ( 2 10 ) ( R nd j R nd j I x R nd R I j R nd j R ref R j a a a σ σ σ σ σ δ σ σ δ δ (6)
For diurnal waves situated far from the resonance, i.e. for ( )

s rad R nd j / 10 . 3 8 - >> -σ σ , we have ) ( R nd j I I j a σ σ δ - ≈
. Therefore, the FCN quality factor is only constrained by the imaginary part of the gravimetric factor of the diurnal waves close to the resonance (K 1 , Ψ 1 and Φ 1 ) and from equation ( 6) we can see that the imaginary part of the resonance strength is mainly constrained by the imaginary part of the gravimetric factors. However, the imaginary parts of the tidal gravity factors are poorly determined, especially for Ψ 1 and Φ 1 that have small amplitudes (particularly a small imaginary part) and are therefore very sensitive to the ocean loading correction error [START_REF] Sato | Parameters of the fluid core resonance inferred from superconducting gravimeter data[END_REF]. We have tested the influence of the error of the imaginary part of the gravimetric factors on the resulting pdfs of the Bayesian inversion. For that, we have considered synthetic gravimetric factors computed from equation ( 2 We increase the relative error on the imaginary part from 1% to 100%, and we also test the case when the error on δ Ι (Ψ 1 ) and δ Ι (Φ 1 ) becomes larger than their values themselves (relative uncertainty larger than 100%). The corresponding pdfs are plotted in Fig. 1. When the relative error on the imaginary part of Φ 1 and Ψ 1 reaches 100%, the pdf for x begins to be non-Gaussian and tends to an asymmetric probability law (Fig. 1 (b)). If we still increase the relative error on δ Ι (Ψ 1 ) and δ Ι (Φ 1 ) to 200% (Fig. 1 (c)) or increase the error for all the waves to 100% (Fig. 1 (e)), then the probability law for x is definitely not Gaussian any more and tends to the infinity towards larger values. For a relative error of 50%, the pdf for x can still be approximated by a Gaussian law (Fig. 1 (d)). In a previous paper, [START_REF] Florsch | Bayesian estimation of the free core nutation parameters from the analysis of precise tidal gravity data[END_REF] obtained a probability law for x that tends to infinity: that is because of the large errors they had on the imaginary parts of the smallest diurnal waves (Ψ 1 and Φ 1 ). As a consequence, if the estimations of the gravimetric factors were more accurate (i.e. with relative errors smaller than 50%), then the probability law for x should be Gaussian. Moreover, using VLBI nutation data, for which the amplitudes of the equivalent Φ 1 and Ψ 1 nutations are high and well-determined, [START_REF] Rosat | FCN resonance parameters from VLBI and Superconducting Gravimeter data[END_REF] have indeed obtained a Gaussian-law for x.

We can also check the influence of the precision on the estimated Ψ 1 gravimetric factor on the FCN resonance parameters. We use the generalized non-linear least-squares formulation of Tarantola and Valette (1982b) and we make the Ψ 1 gravimetric factors varying around its observed value at Strasbourg ( = δ ranges from 0.002 to 0.012. Therefore the uncertainty on the estimation of Q is large even if the ocean loading correction was perfect. The error on Q is therefore mostly due to the fitting error of Ψ 1 because of its small amplitude. There is also a strong correlation between I 1 Ψ δ and a I . Similar computations could be done for K 1 and Φ 1 .

We have shown that the diurnal waves close to the resonance are very sensitive to the FCN parameters. Besides, some correlations exist between the FCN resonance parameters. As

shown by the tilted shapes in Fig. 1, two strong correlations exist between a R and T, on the one hand, and between a I and x, on the other hand. Physically it comes from the response of the whole Earth to the FCN that depends also on the complex eigenfrequency of the core oscillation. Indeed, the damping of the nearly-diurnal free wobble introduced through the complex eigenfrequency could be due to the dissipation by viscomagnetic and topographic coupling at the CMB and to the anelasticity [START_REF] Mathews | Modeling of nutation and precession: New nutation series for nonrigid Earth and insights into the Earth's interior[END_REF] while the imaginary part introduced in the resonance strength, reflects the anelasticity of the Earth which is related to the imaginary part of the Love numbers. As a consequence, the Q-factor is physically correlated to the imaginary part of the resonance strength through the anelasticity of the interface. [START_REF] Florsch | Bayesian estimation of the free core nutation parameters from the analysis of precise tidal gravity data[END_REF] have also shown that the reference gravimetric factor is correlated to a R and hence to the real part of 1 δ . As we have an uncertainty (usually less than 0.3%) on the SG scale factor, this error will be propagated to the estimation of 1 δ . As we will The exploration performed by [START_REF] Florsch | Bayesian estimation of the free core nutation parameters from the analysis of precise tidal gravity data[END_REF] has shown that the imaginary part of the strength I a is close to zero with an error much larger than the value itself. They proposed to cancel this parameter in order to restrict the scanning to a 3-parameter space. However, a I is strongly correlated to Q particularly for Q-values less than 1000: putting a I = 0 forces Q to be around 2600 (cf. Fig. 2 (a)). Note that in this case (elastic model), the obtained T value is close to 465 days, which is the value inferred by [START_REF] Sasao | A simple theory on the dynamical effects of a stratified fluid core upon nutational motion of the Earth[END_REF] for an elastic Earth with a hydrostatic CMB flattening. which is in agreement with predictions of the MHB model [START_REF] Mathews | Modeling of nutation and precession: New nutation series for nonrigid Earth and insights into the Earth's interior[END_REF] and with VLBI nutation observations (e.g. [START_REF] Lambert | The Earth's core parameters as seen by the VLBI[END_REF][START_REF] Koot | Estimation of Earth interior parameters from a Bayesian inversion of very long baseline interferometry nutation time series[END_REF]. In the next part, we will finally estimate the FCN resonance parameters.

To conclude this part, we have seen that the FCN resonance damping is strongly correlated to the anelasticity of the mantle through the imaginary part of the internal Love number. The large errors on the determined phases of the diurnal waves close to the resonance make the parameter x follow a non-Gaussian probability law.

Combined GGP inversion

Here we propose to analyze, in the spectral domain, 7 tidal European SG observations to invert the FCN parameters. Contrary to [START_REF] Ducarme | European tidal gravity observations: Comparison with Earth Tides models and estimation of the Free Core Nutation (FCN) parameters[END_REF], we do not compute averaged gravimetric factors but combine them by using relation (4) and we suppose that the resonance strength a ~ is the same for all the sites. Because a ~corresponds to the mantle response to the The SG sites considered here are Bad-Hambourg (H1, Germany), Moxa (M1, Germany), Membach (MB, Belgium), Medicina (MC, Italy), Strasbourg (ST, France), Vienna (VI, Austria) and Wettzell (W1, Germany). H1 (resp. M1, W1) refer to the lower sphere of the double-sphere SG installed at these German sites. The record length that has been used to retrieve the gravimetric factors from the ETERNA tidal analysis is larger than 5 years for every time-series. These sites have been chosen because they are relatively far from the ocean (except Membach) and their SGs are known to be well calibrated (better than 0.3% accuracy).

The raw data have been corrected for gaps, spikes, steps and other disturbances so that a tidal analysis with ETERNA software [START_REF] Wenzel | The nanogal software: Earth tide data processing package[END_REF] is possible. Before the tidal analysis, the minute data are decimated to 1h (using a filter with a cut-off period of 3h). The ETERNA software then performs a least-squares fit to tides, local air pressure and instrumental drift to retrieve the complex gravimetric factors, the residual gravity, an adjusted barometric admittance, and a polynomial drift function.

The data to be inverted are the complex gravimetric factors corrected for the ocean tide loading effect according to FES2004 ocean model [START_REF] Lyard | Modelling the global ocean tides: modern insights from FES2004[END_REF], for ocean loading computation see for instance [START_REF] Llubes | A multi-technique monitoring of ocean loading in North of France[END_REF]. We have attributed a nominal error of 0.02 nm/s 2 on the ocean loading correction vector. This value has been roughly estimated based on the comparison between different ocean tide models. The inversion is carried out for the four parameters (

I R R nd a a x , , ,σ
) on a 101 4 points space. The explored parameter space is defined based on the results by [START_REF] Florsch | Bayesian estimation of the free core nutation parameters from the analysis of precise tidal gravity data[END_REF], namely x spans from 3 to 8, T varies between 380 and 520 sidereal days, a R is ranging from 0.0004 to 0.001 deg/h and a I varies Based on the MHB Earth's parameter values and using Equation (3), this resonance strength estimate leads to an observed determination of the 1 δ internal pressure gravimetric factor of 1 δ = 0.0306 -i 0.00223 with an error of 2 10 -4 . Theoretical values of the rotational pressure gravimetric factor 1 δ have been computed by [START_REF] Dehant | Analytical approach to the computation of the Earth, the outer core and the inner core rotational motions[END_REF]. A first experimental determination of 1 δ has been proposed by Hinderer et al. (1991b) The problem of negative Q-values encountered in previous studies has been avoided using the logarithm of Q in the inversion scheme.

Conclusion

We have solved the FCN resonance damping equation for four parameters using on the one hand, the more traditional linearized least-squares method optimized with the Levenberg-Marquardt algorithm, and on the other hand, a statistical Bayesian approach. We have demonstrated the good agreement between both results.

Because of the correlation existing between the quality factor Q and the imaginary part of the resonance strength a I , and the necessity for the mantle to possess some anelasticity, this latter parameter cannot be neglected and must be considered in the inversion scheme. Because of the large error on the phase of the diurnal waves close to the resonance, the parameter x = log 10 (Q) does not follow a Gaussian law.

From stacking the SG tidal gravity data, we can improve the determination of the FCN resonance parameters, which are in agreement with VLBI nutation determination. In order to further improve these estimates, the ocean loading correction must be improved at these diurnal frequencies. The nutation has the advantage to be less affected by the ocean loading uncertainty as the contribution of the ocean to the Earth's wobble is much smaller than its effect (mainly local) on the surface gravity.

Time-varying gravity data appear to provide complementary information to the VLBI nutation data as the latter cannot determine the precession while gravity can determine the K 1 tidal wave. Besides, gravity observations enable to retrieve additional Earth's interior parameters like the internal pressure Love numbers.

Table captions

Table 1: A summary of various estimates of period and quality factor of the FCN. In addition

  strength corresponding to the response of the whole Earth to the FCN. The quantity ref δ ~is the value of the gravimetric factor without any resonance process (classical tidal gravimetric factor); it is also the asymptotic value of δj for frequencies far away from the resonance frequency. The eigenperiod T of the FCN expressed in sidereal days in the rotating frame is related to in degrees/solar hour. The quality factor Q,

  e f is the dynamical flattening of the CMB and β the compliance defined by number at the CMB associated with the fluid pressure on the mantle caused by the relative rotation.

  th refers to the theoretical value, and ∆δ is the error on δ (standard deviation).

  based on VLBI results, e.g.[START_REF] Mathews | Modeling of nutation and precession: New nutation series for nonrigid Earth and insights into the Earth's interior[END_REF], a R ≈ 3.4 10 -9 rad/s ≈ 7 10 -4

  ) with x = 4, days, a R = 6 10 -4 deg/h and a I = -5 10 -5 deg/h. Then we assume an error of 0.1% on the real part of the gravimetric factors for nine diurnal waves (from Q 1 to OO 1 waves).

1 Ψδ

 1 attributed as an error on both the real and imaginary parts of 1 ~ψ δ does not take into account the ocean loading correction error as it must reflect the minimal error that we have on the A c c e p t e d M a n u s c r i p t 12 determination of the gravimetric factors at the site. The results of the least-squares inversion lead to the conclusion that, with an error of 0.004 on the real part of 1 ~ψ δ the a posteriori error on T from the least-squares inversion is about 2.3 days. For x = log 10 (Q), the correlation with I is much stronger and x varies between 4.2 and 4.6, corresponding to Q-values between 15849 and 39811, when I 1 Ψ

  next part, this scaling error of 0.3% is smaller than the uncertainty that we obtain for 1 δ .

  When imposing a I equal to a value close to the estimation from the 4-D inversion, i.e. a I = -5 10 -5 deg/h, then we obtain a most probable values for Q of 21241 +/-6956 (cf. Fig. 2 (b))

  a homogeneous response at all SG sites enables to retrieve the global parameters (cf. paragraph 2) like the internal pressure Love number at the CMB. The weight in equation 4 is the combination of the standard deviation of the fitted gravimetric factor by ETERNA software and of the ocean loading error.

  10 -4 deg/h. Therefore we have defined a priori locally uniform laws for these parameters.The pdfs obtained from the Bayesian estimation are plotted in Fig.3. The mean values of the FCN parameters computed from the marginal probability laws are T = 428 +/-3 days, a R = (0.667 +/-0.005) 10 -3 deg/h and a I = (-0.492 +/-0.05) 10 -4 deg/h. The law for Q is not symmetric (and not Gaussian) so we define confidence intervals: Q is between 7762 and 31989 within 90% C.I. Note the good agreement with the Levenberg-Marquardt least-squares results represented by the vertical dotted lines.

  the FCN frequency, the combination of the dynamic flattening of the fluid core and the compliance (or the displacement Love number f h 1

  to theoretical results relative to an elastic Earth, to a slightly inelastic one and to MHB2000 model of[START_REF] Mathews | Modeling of nutation and precession: New nutation series for nonrigid Earth and insights into the Earth's interior[END_REF], we have added experimental results from the IDA (International Digital Accelerometers) network of spring gravimeters and from VLBI (Very Long Baseline Interferometry). The other results are from superconducting gravimeter (SG) datasets: B = Brussel (Belgium), BH = Bad Homburg (Germany), ST = Strasbourg (France), CA = Cantley (Canada), J = 3 Japanese stations, ES = Esashi (Japan), MA = Matsushiro (Japan), CB = Canberra (Australia), MB = Membach (Belgium).

Fig. 1

 1 Figure captions

  and a I used to compute the synthetic delta-factors: x = 4, T = 430 days, a R = 6 10 -4 deg/h and a I = -5 10 -5 deg/h.

Fig. 2

 2 Fig. 2 Joint and marginal pdfs for the FCN parameters (a R , T and x) estimated from Strasbourg SG tidal gravity factors using the Bayesian method with (a) a I = 0 deg/h; (b) a I = -5 10 -5 deg/h. The Marquardt-least squares inversion results are also indicated as vertical

Fig. 3

 3 Fig. 3 Joint and marginal pdfs for the FCN parameters (a R , a I , T and x) estimated from

Fig

  Fig. 1

  1, P 1 , K 1 , Ψ 1 , Φ 1 , J 1 and OO 1 ; (b) 10% for the 7 waves Q 1, O 1 , M 1, P 1 , K 1 , J 1 , OO 1 and 100% for Ψ 1 and Φ 1 ; (c) 10% for Q 1, O 1 , M 1, P 1 , K 1 , J 1 , OO 1 and 200% for Ψ 1 and Φ 1 ; (d) 50% for Q 1, O 1 , M 1, P 1 , K 1 , Ψ 1 , Φ 1 , J 1 and OO 1 ; (e) 100% for Q 1, O 1 , M 1, P 1 , K 1 , Ψ 1 , Φ 1 , J 1 andOO 1 . The vertical dot lines indicate the values of x, T, a R

Acknowledgements

We acknowledge GGP members (http://www.eas.slu.edu/GGP/ggphome.html) for providing SG data. We thank also two anonymous referees for their comments on the manuscript.