Jonas Azzam

Guy David

Tatiana Toro

Wasserstein Distance and the Rectifiability of Doubling Measures: Part II

Keywords: . Rectifiability, Tangent measures, Doubling measures, Wasserstein distance Contents

research documents, whether they are published or not. The documents may come   L'archive ouverte pluridisciplinaire

3 Tangent measures are self-similar 9

1 Introduction

Statement of Results

In this paper we are concerned with understanding the rectifiability properties of doubling measures. Our ultimate goal is to characterize rectifiable doubling measures. Recently Tolsa provided such a characterization for 1-rectifiable measures with upper density bounded below (see [START_REF]Rectifiable measures, square functions involving densities, and Cauchy transform[END_REF]). His conditions are expressed in terms of the properties of the density ratios. We are interested in whether self-similarity properties yield some sort of rectifiabilty. Roughly speaking we analyze how the distance between the dilations of a measure appropriately scaled yield information about the structure of its support. We provide a criteria which ensures that the support of a doubling measure can be decomposed as a union of rectifiable pieces of different dimensions. In a previous paper [START_REF] Azzam | Wasserstein distance and rectifiability of doubling measures: Part I[END_REF], similar decompositions were obtained by looking at conditions that were expressed in terms of the properties of the local distance between the measure and flat measures (that is multiples of Hausdorff measures restricted to affine subsets of Euclidean space). In both cases a minor variant of the L 1 Wasserstein distance is used to estimate the good features of a measure.

To present our results we need to define local distances between measures as well as several quantities which describe the self similar character of a measure. In this paper, µ denotes a Radon measure on R n (i.e., a locally finite positive Borel measure), and Σ = Σ µ denotes its support. That is, Σ = x ∈ R n ; µ(B(x, r)) > 0 for r > 0 ,

where B(x, r) denotes the open ball centered at x and with radius r. We say that µ is doubling when there is a constant C δ > 0 for which µ(B(x, 2r)) ≤ C δ µ(B(x, r)) for all x ∈ Σ and r > 0.

(1.2)

Let B = B(0, 1) denote the unit ball in R n . For M ≥ 0, denote by Lip M (B) the set of functions ψ : R n → R that are M -Lipschitz, i.e., such that

|ψ(x) -ψ(y)| ≤ M |x -y| for x, y ∈ R n , (1.3)
and for which ψ(x) = 0 for x ∈ R n \ B.

(1.4)

Definition 1.1. Let µ and ν be measures on R n , whose restrictions to B := B(0, 1) are probability measures. We set

W 1 (µ, ν) := sup ψ∈Lip 1 (B)
ˆψdµ -ˆψdν .

(1.5)

Thus W 1 (µ, ν) only measures some distance between the restrictions to B of µ and ν. This quantity is similar to the usual L 1 -Wasserstein distance, which by the Kantorovich duality theorem has the same definition as W 1 except that the infimum ranges over all 1-Lipschitz nonnegative functions in B. Note that W 1 has appeared before in the study of rectifiability of measures; see for example [START_REF] Preiss | Geometry of measures in R n : distribution, rectifiability, and densities[END_REF], [START_REF] Tolsa | Uniform rectifiability, Calderón-Zygmund operators with odd kernel, and quasiorthogonality[END_REF], [START_REF]Mass transport and uniform rectifiability[END_REF], and [START_REF] Azzam | Wasserstein distance and rectifiability of doubling measures: Part I[END_REF]. In Section 5, we replace W 1 with a smoother version of local distance W ϕ which is easier to manipulate. Lemmas 5.1 and 5.3 state that W 1 and W ϕ are essentially comparable. We refer to [START_REF] Villani | Optimal transport, Grundlehren der Mathematischen Wissenschaften[END_REF] for a detailed introduction to Wasserstein distances and their properties.

To estimate the self-similarity properties of µ we use several groups of affine transformations of R n . Denote by R the group of affine isometries of R n (i.e., compositions of translations, rotations, and symmetries). Then let G denote the group of similar affine transformations, defined by G = λR ; λ > 0 and R ∈ R .

(1.6)

For G ∈ G , we denote by λ(G) the unique positive number such that G = λ(G)R for some R ∈ R. Denote by D the group of translations and dilations:

D = λI + a ; λ > 0 and a ∈ R n , (1.7)
where I denotes the identity on R n . The transformations that map a given x ∈ R n to the origin, are denoted by To each G ∈ G , we associate the measure µ G = G ♯ µ, which is defined by µ G (A) = µ(G -1 (A)) for every Borel set A ⊂ R n .

G (x) = G ∈ G ; G(x
(1.9)

When G ∈ G (x) for some x ∈ Σ, be may normalize µ G and set

µ G 0 = µ G µ G (B) = µ G µ(G -1 (B)) = µ G µ(B(x, λ(G) -1 ))
(1.10) because µ(B(x, λ(G) -1 )) > 0. This normalization is needed if we want to compute W 1distances.

A special case of this is when G = T x,r , the element of D that maps B(x, r) to B; then µ G and µ G 0 are denoted by µ x,r and µ x,r 0 respectively. That is, µ x,r (A) = µ(x + rA) for A ⊂ R n , and µ x,r 0 = µ x,r µ(B(x, r))

. (1.11) To measure the self-similar nature of µ we introduce two quantities α G and α D . We fix two parameters

1 < λ 1 < λ 2 < ∞. Set G (x, r) = G ∈ G (x) ; λ 1 r ≤ λ(G) -1 ≤ λ 2 r
(1.12) and then α G (x, r) = inf W 1 (µ G 0 , µ x,r 0 ) ; G ∈ G (x, r) .

(1.13) Thus, if α G (x, r) is small, this means that in B, µ x,r 0 is close to some measure µ G 0 , obtained via a transformation G that contracts more than T x,r and possibly rotates as well. After composition with T -1

x,r , the fact that α G (x, r) is small can be interpreted as saying that in B(x, r), µ is quite close to the measure aµ G ′ , where

G ′ = T -1
x,r • G is a contracting element of G that fixes x and a > 0 is a normalizing constant. It is important to note that even though we allow some flexibility in the choice of G and G ′ , we demand that G ′ (x) = x. This is the reason why the usual fractal measures do not satisfy the conditions below.

We also use the analogue of α G (x, r) for the smaller group D.That is, set

D(x, r) = D ∈ D(x) ; λ 1 r ≤ λ(G) -1 ≤ λ 2 r (1.14) and α D (x, r) = inf W 1 (µ G 0 , µ x,r 0 ) ; G ∈ D(x, r) . (1.15)
For α D (x, r) we only compare µ with its image by some dilation centered at x. Obviously α D (x, r) ≥ α G (x, r). Thus conditions on α D (x, r) are more restrictive than those on α G (x, r).

Our goal is to get a control on the part of Σ (see (1.1)) where either α G (x, r), or α D (x, r), are sufficiently small. More precisely, we want to control the sets

Σ 1 = x ∈ Σ ; ˆ1 0 α D (x, r) dr r < ∞ (1.16)
and

Σ 2 = x ∈ Σ ; ˆ1 0 α G (x, r) log(1/r)dr r < ∞ .
(1.17) Theorem 1.2. Let µ be a doubling measure on R n , and denote by Σ its support. Let 1 < λ 1 < λ 2 < ∞ be given, and define the sets Σ 1 and Σ 2 as above. Then there are sets S 0 , ..., S n ⊂ Σ, such that

µ (Σ 1 ∪ Σ 2 ) \ n d=0 S d = 0, (1.18)
and moreover

• S 0 is the set of points where Σ has an atom; it is at most countable, and every point of S 0 is an isolated point of Σ.

• For 1 ≤ d ≤ n, if x ∈ S d
, there exists a d-dimensional vector space V x such that all the tangent measures to µ at x (defined below) are multiples of the Lebesgue measure on V x .

• For 1 ≤ d ≤ n, S d is d-rectifiable, and it can be covered by a countable family of Lipschitz graphs of dimension d.

We may see Theorem 1.2 as a structural decomposition of the good parts of Σ. Tangent measures will play an important role in proof and the definition of the S d . Recall that the set of tangent measures to µ at x, which will be denoted by Tan(µ, x), is the set of non-zero Radon measures σ for which there exist a sequence {r k }, with lim k→∞ r k = 0, and a sequence {a k } of nonnegative numbers, such that σ is the weak limit of the measures a k µ x,r k , (1.19) where the µ x,r k are as in (1.11). That is, for every continuous function f with compact support,

ˆf dσ = lim k→∞ a k ˆf dµ x,r k . (1.20)
Note that since here µ is doubling, Tan(µ, x) is not empty (see for instance the proof of Lemma 2.1 in [START_REF] Azzam | Wasserstein distance and rectifiability of doubling measures: Part I[END_REF]). Furthermore if µ satisfies (1.2) and σ ∈ Tan(µ, x) then σ is also doubling with a constant at most C 2 δ . A priori Tan(µ, x) may be large. Nevertheless Theorem 1.2 ensures that for x ∈ S d , Tan(µ, x) is of dimension 1.

A Lipschitz graph of dimension d is a set Γ A such that

Γ A = x + A(x) ; x ∈ V ,
where V is a vector space of dimension d, A : V → V ⊥ is a Lipschitz map and V ⊥ denotes the (n -d)-dimensional vector space perpendicular to V . In the statement of Theorem 1.2, S d can be covered by Lipschitz graphs where the corresponding function A has Lipschitz constant less than ε, where ε > 0 is any small number given in advance. Note that this yields that S d is d-rectifiable while providing additional information in the sense that S d is completely covered by Lipschitz graphs not simply up to a set of H d -measure zero. Let us make a few more remarks on Theorem 1.2 and its proof. The advantage of using the quantities α G and α D is that they yield information not only about the geometry of the support but also about how the measure is distributed on it. The decomposition of Σ 1 ∪ Σ 2 into pieces of different dimensions is possible once we prove that for µ-almost every point x ∈ Σ 1 ∩ Σ 2 , Tan(µ, x) is entirely composed of flat measures of a same dimension depending on x. Recall that flat measures are multiples of Lebesgue measures on vector subspaces of R n ; that is, for each integer d ∈ [0, n], set

F d = cH d V ; c ≥ 0 and V ∈ G(d, n) , (1.21)
where H d denotes the d-dimensional Hausdorff measure (see [START_REF] Mattila | Geometry of sets and measures in Euclidean spaces[END_REF] or [START_REF] Federer | Geometric measure theory[END_REF]) and

G(d, n) is the set of d-planes in R n . The set of flat measures is F = 0≤d≤n F d .
It is natural to use the self-similarity properties of µ to get information on the structure of Σ, as in Theorem 1.2. In particular the numbers α G (x, r) provide an intrinsic way to measure the regularity of µ. We contrast this approach with the one taken in [START_REF] Azzam | Wasserstein distance and rectifiability of doubling measures: Part I[END_REF] where we were interested on the local approximation of the measure by flat measures.

As we shall see in Section 8, the additional logarithm in (1.17) is used to sum a series which allows us to control the density of µ on most of Σ 2 . It may well be an artifact of the proof.

To prove Theorem 1.2 we need to find a set Σ 0 which covers almost all Σ 1 ∪ Σ 2 and such that all tangents to µ at points in Σ 0 are flat. To accomplish this we define an average analogue of the numbers α G (x, r) by

α * G (x, r) = B(x,r) 2r r α G (y, t)dµ(y)dt (1.22)
(where ffl is our notation for an average) and consider the set

Σ 0 = x ∈ Σ ; lim r→0 α * G (x, r) = 0 . (1.23) Theorem 1.3.
Let µ be a doubling measure on R n , and denote by Σ its support. Let 1 < λ 1 < λ 2 < ∞ be given, and define the functions α G (x, r) and α * G (x, r) and the set Σ 0 as above (see (1.13), (1.22), and (1.23)). Then Tan(µ, x) ⊂ F for every x ∈ Σ 0 .

(1.24)

A consequence of (1.24) and the fact that elements of different F d are far away from each other is that for each x ∈ Σ 0 , there is an integer d ∈ [0, n] such that Tan(µ, x) ⊂ F d .

(1.25) This is not too hard to prove. In the case of Theorem 1.2, we obtain more than (1.25) directly, thus we omit the proof of this fact.

To deduce Theorem 1.2 from Theorem 1.3, we shall first check that

µ (Σ 1 ∪ Σ 2 ) \ Σ 0 = 0. (1.26)
This uses standard techniques from measure theory including the Lebesgue density theorem. Then we show that for each

x ∈ Σ 0 ∩ (Σ 1 ∪ Σ 2 ), Tan(µ, x) = cσ ; c > 0 for some σ ∈ F . (1.27)
Let us now say a few words about the definition of the S d . Set

S d = x ∈ Σ 0 ∩ (Σ 1 ∪ Σ 2 ) ; Tan(µ, x) ⊂ F d ;
(1.28) these sets are disjoint, and by (1.27) or (1.25)

Σ 0 ∩ (Σ 1 ∪ Σ 2 ) = n d=0 S d . (1.29)
The special set S 0 is easily dealt with at the beginning of Section 7, and Theorem 1.2 follows as soon as we prove that for d ≥ 1, S d can be covered by a countable collection (1.30) of Lipschitz graphs of dimension d.

The fact that information on the tangent measures may imply rectifiability properties for the measure is much better understood since D. Preiss [START_REF] Preiss | Geometry of measures in R n : distribution, rectifiability, and densities[END_REF] showed that if µ is a Radon measure, not necessarily doubling, such that for µ-almost every x ∈ R n , the d-density lim r→0 µ(B(x,r)) r d exists and is positive and finite, then R n may be covered, up to a set of µ-measure zero, by a countable collection of d-dimensional Lipschitz graphs. He deduced this from the hypothesis on the d-density and the fact that at µ-almost every x ∈ Σ, Tan(µ, x) ⊂ F . In our case, we are unable to use [START_REF] Preiss | Geometry of measures in R n : distribution, rectifiability, and densities[END_REF] because we are not given any information on the density of µ.

We shall use the fact that since µ is doubling, (1.27) implies the existence of a tangent d-plane to Σ at x, and then (1.30) for the set where (1.27) holds. We include a proof of these simple observations in Section 7.

To prove (1.27), we shall use the numbers α D (x, r) and α G (x, r) to control the variations of the measures µ x,r on Σ 1 and Σ 2 . Eventually we compare them to the tangent measures.

For points of Σ 1 we use the triangle inequality and the summability of the α D (x, r), to show that the distance between two different tangent measures at x is controlled by integrals that tend to 0. To deal with some of the technical complications that arise with the distance W 1 , we shall introduce in Section 5 a smoother variant W ϕ of the Wasserstein distance, study it briefly, and then use it in Section 6 to prove (1.27) on Σ 1 ∩ Σ 0 .

For points of Σ 2 , we'll use the bounds on the numbers α G (x, r) to compute the W ϕdistance between the µ x,r and the tangent measures. This time we can only work modulo rotations, but this is enough to control the W 1 -distance from the µ x,r 0 to flat measures, and apply Theorem 1.5 in [START_REF] Azzam | Wasserstein distance and rectifiability of doubling measures: Part I[END_REF]. This yields additional information on Σ 2 ∩ Σ 0 . In particular, it guarantees that on the sets Σ 2 ∩ S d , µ is absolutely continuous with respect to the Hausdorff measure H d , with a density that can be computed from the measure of balls, and that some local mutual absolute continuity of µ and H d |Σ 2 ∩S d holds. See near (8.17) for a statement, and the rest of Section 8 for the proof.

There is a significant difference between (1.24) (or even (1.25)) and the stronger (1.27). For instance, let Σ be an asymptotically flat snowflake in R 2 , constructed in the usual way but with angles that tend slowly to 0. Put on Σ the natural measure µ, coming from the parameterization of Σ (see [START_REF] David | Asympotically Optimally Doubling Measures and Reifenberg Flat Sets with Vanishing Constant[END_REF]). In this case for µ-almost every x ∈ Σ, Tan(µ, x) = F 1 . Of course Σ is not rectifiable, and Theorem 1.2 says that Σ 1 ∪ Σ 2 is µ-negligible.

The definitions (1.22) and (1.23) ensure that for x ∈ Σ 0 , the local self-similarity character of µ improves as the balls get smaller and smaller, which yields self-similar tangent measures at all point of Σ 0 . That is, we show that if x ∈ Σ 0 , σ ∈ Tan(µ, x) and y lies in the support of σ, there is a transformation G ∈ G such that λ 1 ≤ λ(H) -1 ≤ λ 2 , H(y) = y, and H ♯ σ = cσ for some c > 0. See Lemma 3.1, in Section 3. Once we prove this, showing that σ is flat is mostly a matter of playing with the invariance properties of the support and the measure; see Section 4.

Control of the averages α G

The main goal of this section is to prove (1.26). To this effect define the set

Σ 3 = x ∈ Σ ; ˆ1 0 α G (x, r) dr r < ∞ . (2.1)
Notice that Σ 1 ∪ Σ 2 ⊂ Σ 3 , by (1.16), (1.17), and because α G (x, r) ≤ α D (x, r). Thus (1.26) follows once we prove that

µ(Σ 3 \ Σ 0 ) = 0. (2.2)
For N > 0 large and k ≥ 2, let

Σ 3 (N ) = x ∈ Σ 3 ∩ B(0, N ) ; ˆ1 0 α G (x, r) dr r ≤ N (2.3)
and

ε k = ˆΣ3 (N ) ˆ2-k+2 2 -k α G (y, r) dµ(y)dr r . (2.4) Then k≥2 ε k ≤ 2 ˆΣ3 (N ) ˆ1 0 α G (x, r) dµ(x)dr r ≤ 2N µ(Σ 3 (N )) < ∞. (2.5)
Choose a decreasing sequence {γ k } such that

lim k→∞ γ k = 0 and k≥2 γ -1 k ε k < ∞. (2.6)
For x ∈ Σ 3 (N ), define auxiliary functions α k by

α k (x) = ˆΣ3 (N )∩B(x,2 -k+1 ) ˆ2-k+2 2 -k α G (y, r) dµ(y)dr r . (2.7)
Consider the bad sets

Z k = x ∈ Σ 3 (N ) ; α k (x) ≥ γ k µ(B(x, 2 -k+1 )) . (2.8)
Our goal is to show that Z k is small. Let X ⊂ Z k be a maximal subset whose points lie at distance at least 2 -k+2 from each other. Thus the balls B(x, 2 -k+2 ), x ∈ X, cover Z k , so by (1.2) and (2.8)

µ(Z k ) ≤ x∈X µ(B(x, 2 -k+2 )) ≤ C 2 δ x∈X µ(B(x, 2 -k+1 )) ≤ C 2 δ γ -1 k x∈X α k (x). (2.9) Since the balls B(x, 2 -k+1 ), x ∈ X, are disjoint, x∈X α k (x) ≤ ε k (compare (2.4) and (2.7)); thus µ(Z k ) ≤ C 2 δ γ -1 k ε k .
We are not interested in the precise bound, but merely the fact that k µ(Z k ) < ∞, from which we deduce that if we set

Z * l = k≥l Z k for l ≥ 2, then lim l→∞ µ(Z * l ) = 0. Thus for µ-almost every x ∈ Σ 3 (N ) there is k x ∈ N such that x ∈ Σ 3 (N ) \ Z k for k ≥ k x .
(2.10)

By the Lebesgue differentiation theorem applied to the doubling measure µ, we have that for µ-almost every x ∈ Σ 3 (N )

lim r→0 µ(B(x, r) ∩ Σ \ Σ 3 (N )) µ(B(x, r) ∩ Σ) = 0; (2.11)
see for instance Corollary 2.14 in [START_REF] Mattila | Geometry of sets and measures in Euclidean spaces[END_REF]. Let x ∈ Σ 3 (N ) satisfy (2.10) and (2.11); we want to estimate α * G (x, r) for r small. Choose

k such that 2 -k ≤ r ≤ 2 -k+1 ; then k ≥ k x for r small. Recall from (1.22) that α * G (x, r) = µ(B(x, r)) -1 ˆy∈B(x,r) 2r r α G (y, t)dµ(y)dt.
(2.12)

We decompose the domain of integration above into two parts and estimate each one separately. By (1.2), (2.7), (2.10), and (2.8),

µ(B(x, r)) -1 ˆΣ3 (N )∩B(x,r) 2r r α G (y, t)dµ(y)dt ≤ 4µ(B(x, 2 -k )) -1 ˆΣ3 (N )∩B(x,r) ˆ2-k+2 2 -k α G (y, t) dµ(y)dt t ≤ 4C δ µ(B(x, 2 -k+1 )) -1 α k (x) ≤ 4C δ γ k . (2.13)
This term tends to 0 when r tends to 0, by (2.6).

For the second part, we notice that α G (y, t) ≤ 2 by definition (see (1.13) and (1.5)), so

µ(B(x, r)) -1 ˆΣ∩B(x,r)\Σ 3 (N ) 2r r α G (y, t)dµ(y)dt ≤ 2µ(B(x, r)) -1 µ(Σ ∩ B(x, r) \ Σ 3 (N )), (2.14)
which tends to 0 by (2.11). Combining (2.13) and (2.14) we get that

lim r→0 α * G (x, r) = 0 for µ-almost every x ∈ Σ 3 (N ).
In other words,

µ(Σ 3 (N ) \ Σ 0 ) = 0 (see (1.23)); (2.
2) follows easily, and so does (1.26).

3 Tangent measures are self-similar

In this section we start the proof of Theorem 1.3 and prove the basic self-similarity estimate for tangent measures.

Lemma 3.1. Let µ be a doubling measure, let Σ denote its support, let Σ 0 ⊂ Σ be as in (1.23), and σ ∈ Tan(µ, x) be a tangent measure of µ at a point x ∈ Σ 0 . For each y ∈ Ξ, the support of σ, there exist H ∈ G such that H(y) = y,

λ 1 ≤ λ(H) -1 ≤ λ 2 (3.1)
and

H ♯ σ = cσ for some c > 0. (3.2)
The numbers λ 1 and λ 2 are the same as in (1.12), the dilation number λ(H) is defined below (1.6), and H ♯ σ, the push forward image of σ by H, is defined as in (1.9).

Proof. We may assume, without loss of generality, that x = 0. Since σ ∈ Tan(µ, x) there are coefficients a k ≥ 0 and radii r k > 0, such that lim k→∞ r k = 0, and σ is the weak limit of the measures {σ k }, where

σ k = a k µ 0,r k = a k µ R k with R k (u) = r -1 k u for u ∈ R n . (3.3) Note that R k maps B(x, r k ) = B(0, r k ) to B. See (1.11) for the definition of µ x,r . Let α k = sup α * G (0, r) ; 0 < r < √ r k ; (3.4) then since x ∈ Σ 0 (see (1.23)) lim k→∞ α k = 0. (3.5)
By (3.4), if for k large r k < ρ k < √ r k (3.6)
then α * G (0, ρ k ) ≤ α k for these k. Since σ is the weak limit of the σ k , for each y ∈ Ξ := supp σ we can find points

y k ∈ supp(σ k ) = r -1 k Σ such that lim k→∞ |y k -y| = 0. (3.7)
Let {η k } k≥1 and {ρ k } k≥1 be sequences such that (3.6) holds for k large, and also

lim k→∞ ρ k r k = ∞ and lim k→∞ η k = 0. (3.8)
Consider

A k = B(r k y k ,η k r k ) 2ρ k ρ k α G (z, t)dµ(z)dt. (3.9)
For k large, the domain of integration B(r

k y k , η k r k ) is contained in B(0, ρ k ) (because y k tends to y, η k tends to 0, and r -1 k ρ k tends to +∞; see (3.7) and (3.8)). Recall that B(0, ρ k )×[ρ k , 2ρ k ] is the domain of integration in the definition of α * G (0, ρ k ) (see (1.22)); hence for k large A k ≤ µ(B(0, ρ k )) µ(B(r k y k , η k r k )) α * G (0, ρ k ) ≤ C 2+log 2 (ρ k /(η k r k )) δ α k (3.10)
by (3.9), (1.22), (3.6), (3.4), and the doubling property (1.2). Later on, we will choose ρ k and η k , depending on α k , so that A k is still small enough.

By Chebyshev's inequality there exist

z k ∈ Σ ∩ B(r k y k , η k r k ) and t k ∈ [ρ k , 2ρ k ] (3.11) such that α G (z k , t k ) ≤ A k . (3.12)
By the definition of α G (see (1.13)) there exists

G k ∈ G (z k , t k ) such that W 1 (µ G k 0 , µ z k ,t k 0 ) ≤ 2A k , (3.13)
which by (1.5) means that

ˆψdµ G k 0 -ˆψµ z k ,t k 0 ≤ 2A k for any ψ ∈ Lip 1 (B). (3.14) Our next goal is to interpret (3.14) in terms of σ k . Let T k ∈ D be such that for u ∈ R n T k (u) = u -z k t k and so T k (B(z k , t k )) = B. (3.15)
The definitions (1.10) and (1.11) yield

µ G k 0 = e k µ G k and µ z k ,t k 0 = e ′ k µ z k ,t k = e ′ k µ T k , (3.16)
where e k and e ′ k come from the normalization, and are given by

e k = µ G k (B) -1 and e ′ k = µ T k (B) -1 . (3.17)
Let ψ be any Lipschitz function supported on B and set I k = e k ´Gk (Σ) ψdµ G k ; by (3.16) and (1.9),

I k = e k ˆGk (Σ) ψdµ G k = e k ˆΣ ψ(G k (ξ))dµ(ξ) = e k ˆψ • G k dµ. (3.18) By (3.3), σ k = a k µ R k = a k (R k ) ♯ µ, hence µ = a -1 k (R -1 k ) ♯ σ k .
Thus a similar computation to the one in (3.18) yields

I k = e k a -1 k ˆψ • G k • R -1 k dσ k . (3.19)
A similar computation, with G k replaced by T k , yields

I ′ k := ˆψµ z k ,t k 0 = e ′ k a -1 k ˆψ • T k • R -1 k dσ k . (3.20)
We want to apply (3.19) and (3.20) to special functions ψ. Let ϕ be a compactly supported 1-Lipschitz function. Let

ψ = ϕ • R k • T -1 k . (3.21)
Note that ψ is a Lipschitz function with constant less or equal to r -1

k t k ≤ 2r -1 k ρ k (see (3.11)). If ϕ is supported in B(0, R), for k large enough the support of ψ is contained in T k • R -1 k (B(0, R)) = T k (B(0, r k R)) = B -z k t k , t -1 k r k R ⊂ B(0, r k t -1 k (R + η k + |y k |)) ⊂ B(0, 1),
where we have used (3.11), (3.8), and (3.7). Because of (3.10), ψ = (2r -1 k ρ k ) -1 ψ is 1-Lipschitz; then (3.14) applies to ψ, and (3.10) yields

|I ′ k -I k | ≤ 4A k r -1 k ρ k ≤ 4r -1 k ρ k C 2+log 2 (ρ k /(η k r k )) δ α k =: α k (3.22)
where I k and I ′ k are as in (3.18) and (3.20) with ψ coming from (3.21). The final identity is the definition of α k . Notice that even though ϕ does not depend on k, ψ does, but this is not an issue.

Note that by (3.21) and (3.20), we have

I ′ k = e ′ k a -1 k ˆϕdσ k . (3.23)
Similarly, by (3.19) and (3.21) we have

I k = e k a -1 k ˆϕ • R k • T -1 k • G k • R -1 k dσ k . (3.24) Set H k = R k • T -1 k • G k • R -1 k . (3.25)
Then H k ∈ G , and by (3.15), its dilation factor λ(H k ) (which is also the n-th root of its Jacobian determinant) is such that

λ(H k ) -1 = λ(T k )λ(G k ) -1 = t -1 k λ(G k ) -1 ∈ [λ 1 , λ 2 ], (3.26)
because

G k ∈ G (z k , t k ),
and by the definition (1.12). Note that by (3.3), the fact that

G k ∈ G (z k , t k ) (see (1.
12) and (1.8)), and (3.15), we have

H k (r -1 k z k ) = R k • T -1 k • G k • R -1 k (r -1 k z k ) = R k • T -1 k • G k (z k ) = R k • T -1 k (0) = R k (z k ) = r -1 k z k . (3.27)
Notice also that by (3.11)

|r -1 k z k -y| ≤ |r -1 k z k -y k | + |y k -y| = r -1 k |z k -r k y k | + |y k -y| ≤ η k + |y k -y|.
e k a -1 k ˆϕ • H k dσ k -e ′ k a -1 k ˆϕdσ k ≤ α k . (3.30)
By (3.17), (1.9), and (3.15),

e ′ k = µ T k (B) -1 = µ(T -1 k (B)) -1 = µ(B(z k , t k )) -1 . (3.31) Similarly notice that G k (z k ) = 0 because G k ∈ G (z k , t k ) ⊂ G (z k ) (see (1.
12) and (1.8)); then (3.17) and (1.9) yield

e k = µ G k (B) -1 = µ(G -1 k (B)) -1 = µ(B(z k , λ(G k ) -1 )) -1 . (3.32)
In addition,

B(z k , λ 1 t k ) ⊂ G -1 k (B) = B(z k , λ(G k ) -1 ) ⊂ B(z k , λ 2 t k ) (3.33) because G k (z k ) = 0 and
a k µ(B(0, r k )) = σ k (B) ≤ ˆf dσ k ≤ σ k (2B) = a k µ(B(0, 2r k )) ≤ C δ a k µ(B(0, r k )).
(3.35) Moreover, since σ is also doubling (see the remark below (1.20)), we have that

σ(B) ≤ ˆf dσ ≤ σ(2B) ≤ C 2 δ σ(B). (3.36)
Thus by (3.35), (3.36) and the definition of σ there exists C > 1 such that for k large,

C -1 a k µ(B(0, r k )) ≤ σ(B) ≤ CC 2 δ a k µ(B(0, r k )). ( 3
(e ′ k ) -1 = µ(B(z k , t k )) ≤ µ(B(0, 4ρ k )) ≤ C 3+log 2 (ρ k /r k ) δ µ(B(0, r k )). (3.38)
Combining (3.37) and (3.38) we obtain

a k ≤ CC 5+log 2 (ρ k /r k ) δ e ′ k • σ(B). ( 3
b k ˆϕ • H k dσ k -ˆϕdσ k ≤ a k α k e ′ k ≤ CC 5+log 2 (ρ k /r k ) δ α k σ(B). (3.40)
Now we choose ρ k and η k . Denote by Since ϕ is Lipschitz and compactly supported, so is ϕ • H and

β k = CC 5+log 2 (ρ k /r k ) δ α k = CC 5+log 2 (ρ k /r k ) δ • 4r -1 k ρ k C 2+log 2 (ρ k /(η k r k )) δ
lim k→∞ ˆϕ • H dσ k = ˆϕ • H dσ (3.43)
because σ is the weak limit of the σ k . Note that there is also a ball B such that for

k large ϕ • H(x) = ϕ • H k (x) = 0 for x ∈ R n \ B; then ˆ|ϕ • H -ϕ • H k | dσ k ≤ ||ϕ|| lip ˆB |H -H k | dσ k ≤ ||ϕ|| lip ||H -H k || L ∞ (B) σ k (B) ≤ 2||ϕ|| lip ||H -H k || L ∞ (B) σ(2B). (3.44)
Thus

lim k→∞ ˆϕ • H k dσ k -ˆϕ • H dσ k = 0. ( 3

Self-similar measures are flat

In this section we complete the proof of Theorem 1.3. Using the notation in Section 3, our goal is to show that if σ ∈ Tan(µ, x 0 ), where x 0 ∈ Σ 0 (see (1.23)) then σ is a flat measure. Lemma 3.1 guarantees that for each y ∈ Ξ (the support of σ), there is a transformation H(y) ∈ G and a constant c(y) > 0 such that

λ(y) := λ(H(y)) ∈ [λ -1 2 , λ -1 1 ] (4.1)
and

H(y) ♯ σ = c(y)σ. (4.2)
By definition of G (see (1.6)), the linear part of H(y) is of the form λ(y)R(y), where R(y) is a linear isometry. Since H(y) fixes y, this means that

H(y)(u) = y + λ(y)R(y)(u -y) for u ∈ R n . (4.
3)

The next lemma allows us to replace H(y) with one of its large powers, chosen so that its isometric part is very close to the identity.

Lemma 4.1. For each choice of ε > 0, there is an integer m 0 , that depends only on ε and n, such that for each y ∈ Ξ and each integer ℓ ≥ 1, there is an integer m(y)

∈ [1, m 0 ] such that ||R(y) m(y)ℓ -I|| ≤ ε. (4.4)
Proof. Here I is the identity mapping. We use the compactness of the group of linear isometries of R n to choose m 0 large enough so that if R 1 , . . . , R m 0 are linear isometries, we can find integers

m 1 , m 2 such that 1 ≤ m 1 < m 2 < m 0 and ||R m 2 -R m 1 || ≤ ε.
We apply this with R m = R(y) mℓ to find m 1 and m 2 such that ||R(y)

m 2 ℓ -R(y) m 1 ℓ || ≤ ε. Then ||R(y) (m 2 -m 1 )ℓ -I|| ≤ ε, as needed.
Next we study elementary properties of Ξ. Notice that by (4.2), H(y)(Ξ) = Ξ, and iterations also yield

H(y) m (Ξ) = Ξ for m ≥ 1. (4.5)
Lemma 4.2. The set Ξ is convex.

Proof. Let x, y ∈ Ξ be given. Our goal is to show that the segment [x, y] is contained in Ξ.

For each ε > 0 and ℓ ≥ 1, we construct a sequence {y k } in Ξ (depending on ε and ℓ) which will allow us to estimate how far [x, y] is from Ξ. We start with y 0 = y. If k ≥ 0 and y k ∈ Ξ has been defined, we define y k+1 as follows. Set

H k = H(y k ) m(y k )ℓ and y k+1 = H k (x). (4.6)
By (4.5) and since x ∈ Ξ, y k+1 ∈ Ξ. Let us show that for ℓ large and ε small, the y k stay close to the segment [x, y] and converge (slowly) to x. First observe that by iterations of (4.3), H k is given by

H k (u) = y k + λ ′ k R k (u -y k ) for u ∈ R n , (4.7)
with

λ ′ k = λ(H(y k )) m(y k )ℓ ∈ [λ -m 0 ℓ 2 , λ -ℓ 1 ] and hence λ ′ k < 1. (4.8)
In addition, R k = R(y k ) m(y k )ℓ and therefore, by (4.4),

||R k -I|| ≤ ε. Set r k = |x -y k | and y * k+1 = y k + λ ′ k (x -y k ).
Notice that by (4.7) and (4.4),

|y k+1 -y * k+1 | = [y k + λ ′ k R k (x -y k )] -[y k + λ ′ k (x -y k )] = λ ′ k [R k -I](x -y k ) ≤ ελ ′ k r k . (4.9) Since |y * k+1 -x| = (1 -λ ′ k )|y k -x| = (1 -λ ′ k )r k , (4.10)
we also get that if ε < 1/2, Using the definition of y * k+1 , (4.9), and (4.8) we can estimate their successive distances Proof. Let V be the smallest affine subspace of R n that contains Ξ, and let d denote its dimension. Choose d + 1 affinely independent points y 0 , . . . y d in V (this means that the vectors y j -y 0 , j ≥ 1, are linearly independent), and set y = 1 d+1 d j=0 y j . By Lemma 4.2, Ξ contains the convex hull of the y j , so there is a small radius r > 0 such that V ∩ B(y, r) ⊂ Ξ.

r k+1 = |y k+1 -x| ≤ |y * k+1 -x| + |y k+1 -y * k+1 | ≤ (1 -λ ′ k )r k + ελ ′ k r k ≤ (1 -λ ′ k /2)r k . ( 4
|y k+1 -y k | ≤ |y k+1 -y * k+1 | + |y * k+1 -y k | ≤ ελ ′ k r k + λ ′ k r k ≤ 2λ ′ k r k ≤ 2λ -ℓ 1 |x -y|. ( 4
Recall from (4.5) that H(y) m (Ξ) = Ξ for m ≥ 1. Applying the bijection H(y) -m to both sides, we see that H(y) -m (Ξ) = Ξ for m ≥ 1, so

H(y) -m (V ∩ B(y, r)) ⊂ Ξ. (4.16)
We know that H(y) -m (V ∩B(y, r)) is a nontrivial open subset of the affine space H(y) -m (V ), and since Ξ ⊂ V we get that H(y) -m (V ) ⊂ V . Then by a dimension count H(y) -m (V ) = V , and by the description (4.3) of H(y), we see that

H(y) -m (V ∩ B(y, r)) = V ∩ B(y, λ(H(y)) -m r). (4.17)
Recall that λ(H(y)) < 1; then λ(H(y)) -m r is as large as we want by picking m as large as we need. Hence (4.16) guarantees that Ξ contains V , and this means that Ξ = V . Note that Remark 14.4 (2) in [START_REF] Mattila | Geometry of sets and measures in Euclidean spaces[END_REF] ensures that 0 ∈ Ξ, as Ξ is the support of σ and σ ∈ Tan(µ, x 0 ) where µ is a doubling measure and x 0 ∈ Σ. Hence Ξ = V is a vector space. This completes our proof of Lemma 4.3. Now we study the distribution of σ on Ξ. If Ξ is reduced to the origin, then σ is a Dirac mass, and Dirac masses lie in F 0 . Thus in this case σ is trivially flat. We may now assume that Ξ is a vector space of dimension d > 0. Proof. Since σ is Ahlfors regular of dimension D (by Lemma 4.4), a standard covering argument (see for instance Lemma 18.11 in [START_REF] David | Singular sets of minimizers for the Mumford-Shah functional[END_REF]) guarantees that there exists a constant C > 0 such that This completes the proof of Theorem 1.3. In fact we have proved that if σ ∈ Tan(µ, x 0 ) with x 0 ∈ Σ 0 (see (1.23)) then σ = c 0 H d Ξ for some vector space Ξ.

C -1 H D Ξ ≤ σ ≤ CH D Ξ, ( 4

A smoother version of the Wasserstein W 1 distance

So far we managed to work with the distance W 1 defined by (1.5), but for the proof of Theorem 1.2, it is more convenient to use a slightly smoother variant, which attenuates the possible discontinuities in r > 0 of the normalizing factors µ x,r (B(0, 1)

) -1 = µ(B(x, r)) -1 .
Let us choose a smooth radial function ϕ such that Recall that Lip 1 (B) is defined near (1.3). The distance W 1 above essentially corresponds to ϕ = 1 B here.

1 B(0,1/2) ≤ ϕ ≤ 1 B(0,
We required (5.2) (and (5.1)) to make sure that we do not divide by 0. But notice that even when µ(B(0, 1/2)) or ν(B(0, 1/2)) is very small, we always get that

W ϕ (µ, ν) ≤ 2,
(5.4) because | ´ψϕdµ| ≤ ´ϕdµ and similarly for ν. Note that W ϕ (aµ, bν) = W ϕ (µ, ν) for a, b > 0, (5.5) so we do not need to normalize µ and ν in advance. Finally observe that W ϕ satisfies the triangle inequality. That is, if σ is a third Radon measure such that σ(B(0, 1/2)) > 0, if follows at once from the definition that

W ϕ (µ, σ) ≤ W ϕ (µ, ν) + W ϕ (ν, σ). (5.6)
Let us check that if we restrict to measures that are not to small on B(0, 1/2), then W 1 controls W ϕ .

Lemma 5.1. Let µ and ν be Radon measures such that (5.2) holds and

µ(B) = ν(B) = 1.
(5.7)

Then W ϕ (µ, ν) ≤ 1 + 2||ϕ|| lip µ(B(0, 1/2)) W 1 (µ, ν), (5.8)
where ||ϕ|| lip denotes the Lipschitz norm of ϕ.

The fact that the estimate is not symmetric is not an issue. In particular we shall apply (5.8) to doubling measures µ and ν; in this case µ(B(0,

1/2)) ∼ µ(B) = 1 = ν(B) ∼ ν(B(0, 1/2)).
Proof. Let ψ ∈ Lip 1 (B) be given. The definition (1.5), applied to ψϕ, yields

ˆψϕdµ -ˆψϕdν ≤ ||ψϕ|| lip W 1 (µ, ν) ≤ (1 + ||ϕ|| lip )W 1 (µ, ν).
(5.9)

The same definition, applied to ϕ itself, yields

ˆϕdµ -ˆϕdν ≤ ||ϕ|| lip W 1 (µ, ν). ( 5
∆ ≤ |a -c| + |d -b| b ≤ 1 + 2||ϕ|| lip µ(B(0, 1/2)) W 1 (µ, ν).
(5.13)

Taking the supremum over ψ ∈ Lip 1 (B) in (5.13) ( recall (5.11)) yields (5.8).

The following lemma specifies the sense in which W ϕ is more stable that W 1 .

Lemma 5.2. Let µ and ν be Radon measures and let θ ∈ (0, 1/2] be such that µ(B(0, θ/2)) > 0 and ν(B(0, θ/2)) > 0, (5.14) and define new measures µ 1 and ν 1 by

µ 1 (A) = µ(θA) and ν 1 (A) = ν(θA) for A ⊂ R n .
(5.15)

Then W ϕ (µ 1 , ν 1 ) ≤ θ -1 (1 + 4||ϕ|| lip ) µ(B(0, 1)) µ(B(0, θ/2))
W ϕ (µ, ν).

(5.16)

As in (5.8) the estimate is not symmetric in µ and ν, but is nonetheless valid. We require that µ(B(0, θ/2)) = 0 and ν(B(0, θ/2)) = 0 to make sure that W ϕ (µ 1 , ν 1 ) is easily defined. Often µ is the restriction to B of a doubling measure and its support contains the origin; then µ(B(0,1)) µ(B(0,θ/2)) ≥ C -1 , for some C that depends on θ and the doubling constant C δ .

Proof. Let ψ ∈ Lip 1 (B) be given; we want to control the quantity

∆ = ´ψϕdµ 1 ´ϕdµ 1 - ´ψϕdν 1 ´ϕdν 1 =: a b - c d (5.17)
(as above, but with integrals relative to µ 1 and ν 1 ). Notice that by (5.15),

a = ˆψϕdµ 1 = ˆψ(θ -1 x)ϕ(θ -1 x)dµ(x) = ˆψ(θ -1 x)ϕ(θ -1 x)ϕ(x) 2 dµ(x), (5.18)
where we just use the fact that ϕ(θ -1 x) = 0 when x ∈ R n \ B(0, 1/2), and the special shape of ϕ in (5.1), to add an extra ϕ 2 (x). Similarly,

c = ˆψϕdν 1 = ˆψ(θ -1 x)ϕ(θ -1 x)ϕ 2 (x)dν(x). (5.19)
The same computations without ψ yield

b = ˆϕdµ 1 = ˆϕ(θ -1 x)ϕ(x) 2 dµ(x), (5.20)
d = ˆϕdν 1 = ˆϕ(θ -1 x)ϕ(x) 2 dν(x). (5.21)
It is also useful to introduce e = ˆϕdµ and f = ˆϕdν.

(5.22)

Let us first estimate δ 1 = a e -c f . We want to apply the definition of W ϕ (µ, ν) to the function Ψ defined by Ψ(x) = ψ(θ -1 x)ϕ(θ -1 x)ϕ(x). Notice that Ψ is supported in B (this is why we added ϕ(x)), and its Lipschitz norm is at most θ

-1 (1 + ||ϕ|| lip ) + ||ϕ|| lip ≤ θ -1 (1 + 2||ϕ|| lip ). Thus (5.3) yields |δ 1 | = a e - c f ≤ θ -1 ||Ψ|| lip W ϕ (µ, ν) ≤ θ -1 (1 + 2||ϕ|| lip )W ϕ (µ, ν). (5.23)
We can also apply the definition of W ϕ (µ, ν) to ϕ(θ -1 x)ϕ(x), whose Lipschitz norm is at most 2θ -1 ||ϕ|| lip , and we get that

|δ 2 | := b e - d f ≤ 2θ -1 ||ϕ|| lip W ϕ (µ, ν). (5.24)
Thus

a b = e b a e = e b c f + δ 1 = e b c d d f + eδ 1 b = e b c d b e + δ 2 + eδ 1 b = c d + ecδ 2 bd + eδ 1 b ,
where δ 1 and δ 2 are as in(5.23) and (5.24). Thus

∆ = a b - c d ≤ |ecδ 2 | bd + |eδ 1 | b . (5.25)
Now |c| ≤ d because ϕ ≥ 0 and |ψ| ≤ 1, e = ´ϕdµ ≤ µ(B(0, 1)), and b = ´ϕ(θ -1 x)ϕ(x) 2 dµ(x) ≥ µ(B(0, θ/2)) because of (5.1). Thus we have

∆ ≤ (δ 1 + δ 2 ) e b ≤ θ -1 (1 + 4||ϕ|| lip )W ϕ (µ, ν) µ(B(0, 1)) µ(B(0, θ/2)) . (5.26)
Noting (5.17) and taking the supremum over ψ ∈ Lip 1 (B), we obtain (5.16).

The next lemma is used in Section 8. It shows that the distance function W ϕ also controls W 1 in some averaged way. Thus W 1 and W ϕ are basically interchangeable. Lemma 5.3. Let µ and ν are Radon measures such that µ(B(0, 1/4)) > 0 and ν(B(0, 1/4)) > 0. Define µ t and ν t , 1/4 ≤ t ≤ 1/2, by

µ t (A) = µ(tA) µ(B(0, t)) and ν t (A) = µ(tA) ν(B(0, t))
for A ⊂ R n .

(5.27)

Then 1/2 1/4 W 1 (µ t , ν t )dt ≤ (8 + ||ϕ|| lip )µ(B) µ(B(0, 1/4)) W ϕ (µ, ν).
(5.28)

Proof. First notice that the statement does not change when we multiply µ and ν by positive constants. So we may assume that ˆϕdµ = ˆϕdν = 1.

(5.29)

Next fix t ∈ [1/4, 1/2] and ψ ∈ Lip 1 (B). We want to estimate

ˆψdµ t -ˆψdν t = ´ψ(t -1 x)dµ(x) µ(B(0, t)) - ´ψ(t -1 x)dν(x) ν(B(0, t)) =: a b - c d .
(5.30)

As before,

∆ = |ad -bc| bd = |d(a -c) -c(b -d)| bd ≤ |a -c| + |b -d| b (5.31) because |c| = | ´ψ(t -1 x)dν(x)| ≤ ||ψ|| ∞ ν(B(0, t)) ≤ d. Notice that ϕ(x) = 1 when ψ(t -1 x) = 0 since this implies that |t -1 x| ≤ 1 and hence |x| ≤ t ≤ 1/2. Thus |a -c| = ˆψ(t -1 x)dµ(x) -ˆψ(t -1 x)dν(x) (5.32) = ˆψ(t -1 x)ϕ 2 (x)dµ(x) -ˆψ(t -1 x)ϕ 2 (x)dν(x) .
We apply the definition (5.3) of W ϕ with the function x → ψ(t -1 x)ϕ(x), which is supported in B and (t -1 + ||ϕ|| lip )-Lipschitz. We obtain that

|a -c| ≤ (t -1 + ||ϕ|| lip )W ϕ (µ, ν) (5.33)
Notice also that

1 b = ´ϕdµ µ(B(0, t)) ≤ µ(B) µ(B(0, t)) ≤ µ(B) µ(B(0, 1/4)) . (5.34)
Thus

∆ ≤ |a -c| + |b -d| b ≤ µ(B) µ(B(0, 1/4)) [(t -1 + ||ϕ|| lip )W ϕ (µ, ν) + |b -d|] = µ(B) µ(B(0, 1/4)) (t -1 + ||ϕ|| lip )W ϕ (µ, ν) + |µ(B(0, t)) -ν(B(0, t))| . (5.35)
We take the supremum over ψ ∈ Lip 1 (B) and get that

W 1 (µ t , ν t ) ≤ µ(B) µ(B(0, 1/4)) (t -1 + ||ϕ|| lip )W ϕ (µ, ν) + |µ(B(0, t)) -ν(B(0, t))| . (5.36) Since t -1 ∈ [2, 4],
(5.28) will follow as soon as we prove that

ˆ[1/4,1/2] |µ(B(0, t)) -ν(B(0, t))|dt ≤ W ϕ (µ, ν). (5.37)
Let h be a bounded measurable function, defined on [1/4, 1/2]; we want to evaluate

I h = ˆ[1/4,1/2] h(t)[µ(B(0, t)) -ν(B(0, t))]dt. (5.38)
Observe that by Fubini

ˆ[1/4,1/2] h(t)µ(B(0, t))dt = ˆx∈B(0,1/2) ˆ1t∈[1/4,1/2] 1 t>|x| h(t)dt dµ(x),
and similarly for ν. Set ψ h (x) = ´1t∈[1/4,1/2] 1 t>|x| h(t)dt. This is a ||h|| ∞ -Lipschitz function of |x|, which vanishes when |x| ≥ 1/2, so by (5.1), (5.3) and the normalization in (5.29) we have

|I h | = ˆB(0,1/2) ψ h dµ - ˆB(0,1/2) ψ h dν (5.39) = ˆψh ϕdµ -ˆψh ϕdν ≤ ||h|| ∞ W ϕ (µ, ν).
Thus (5.39) holds for all bounded measurable functions h defined on [1/4, 1/2], and (5.37) follows by duality. We saw earlier that (5.37) implies (5.28). Lemma 5.3 follows.

We conclude this section with an easy observation concerning the behavior of W ϕ (µ, ν) when taking limits.

Lemma 5.4. Let µ and ν satisfy (5.2), suppose that µ is the weak limit of some sequence {µ k }, and that ν is the weak limit of some sequence {ν k }. Then

W ϕ (µ, ν) ≤ lim inf k→∞ W ϕ (µ k , ν k ).
(5.40)

Proof. Set L k = W ϕ (µ k , ν k ) and L = lim inf k→∞ L k .
Notice that ´ϕdµ = lim k→∞ ´ϕdµ k by weak convergence; by (5.2), this implies that ´ϕdµ k > 0 for k large. The same argument applied to a continuous function f ≤ 1 B(0,1/2) such that ´f dµ > 0 shows that µ k (B(0, 1/2)) > 0 for k large. Similar observations hold for ν and ν k . For each ψ ∈ Lip 1 (B), the weak convergence yields ´ψϕdµ = lim k→∞ ´ψϕdµ k . For k large,

´ψϕdµ k ´ϕdµ k - ´ψϕdν k ´ϕdν k ≤ L k .
Since each term has a limit and the denominators are bounded away from 0, taking a lim inf we have that ´ψϕdµ ´ϕdµ -´ψϕdν ´ϕdν ≤ L.

Taking the supremum over ψ ∈ Lip 1 (B) we conclude that (5.40) holds.

6 Uniqueness of the tangent measure at "good points"

In this section show that for x ∈ Σ 0 ∩ Σ 1 , Tan(µ, x) is a one-dimensional set of flat measures.

Recall that Σ 1 and Σ 0 were defined in (1.16) and (1.23) respectively. Lemma 6.1. Let µ be a doubling measure, and let x ∈ Σ 0 ∩ Σ 1 . Then there is a nonzero flat measure σ such that Tan(µ, x) = cσ ; c > 0 .

Proof. Fix µ and x ∈ Σ 0 ∩ Σ 1 ; without loss of generality, we may assume that x = 0. By Theorem 1.3, we know that Tan(µ, 0) ⊂ F (6.1)

where F denotes the set of flat measures (see (1.21)). By definition of Σ 1 ,

ˆ1 0 α D (0, r) dr r < ∞. (6.2)
Thus it only remains to show that Tan(µ, 0) is one-dimensional. Our initial goal is to bound the W ϕ distance for two different scaled dilations of µ by α D (0, •) at the right scale. For each r ∈ (0, 1/4) we use Chebyshev's inequality to find r + ∈ [2r, 4r] such that

α D (0, r + ) ≤ (2r) -1 ˆ4r 2r α D (0, t)dt ≤ 2 ˆ4r 2r α D (0, t) dt t . (6.3)
By the definition of α D (0, r + ) (see (1.15)): there is a transformation G ∈ D(0, r + ), such that

W 1 (µ G 0 , µ 0,r + 0 ) ≤ 2α D (0, r + ). (6.4)
By the definition (1.14) of D(0, r + ), G is simply given by G(u) = λu, with λ

-1 = λ(G) -1 ∈ [λ 1 r + , λ 2 r + ]. Set r * = λ -1 ∈ [λ 1 r + , λ 2 r + ]; (6.5)
notice that G is the homotety that sends B(0, r * ) to B, so µ G 0 = µ 0,r * 0 (see near (1.11)) and now (6.4) says that W 1 (µ 0,r * 0 , µ

0,r + 0 ) ≤ 2α D (0, r + ). (6.6)
First apply Lemma 5.1 to the measures µ 0,r * 0 and µ 0,r + 0

;(6.6) yields that W ϕ (µ 0,r * 0 , µ 0,r + 0

) ≤ Cα D (0, r + ) (6.7) (where we do not record any more the dependence on ϕ or C δ ). Then we apply Lemma 5.2 to the same measures, with θ = r/r * . Notice that θ ≤ 1/2 because r * ≥ r + ≥ 2r.

Recall from (5.15) that µ 1 is defined by

µ 1 (A) = µ 0,r * 0 (θA) = µ 0,r * 0 (rA/r * ) = µ(r * rA/r * ) µ(B(0, r * )) = µ(rA) µ(B(0, r * )) = µ(B(0, r)) µ(B(0, r * )) µ 0,r 0 (A) (6.8)
by (1.11). Similarly, ν 1 is defined by

ν 1 (A) = µ 0,r + 0 (θA) = µ 0,r + 0 (rA/r * ) = µ(r + rA/r * ) µ(B(0, r + )) = µ(ρ(r)A) µ(B(0, r + )) = µ(B(0, ρ(r))) µ(B(0, r + )) µ 0,ρ(r) 0 (A) (6.9) where ρ(r) = rr + r * ∈ [λ -1 2 r, λ -1 1 r]. (6.10)
By (5.5), (6.8), (6.9), Lemma 5.2, the fact that µ 0,r * 0 is also doubling with the same constant as µ (which controls mass ratio in (5.16)), and (6.7), W ϕ (µ 0,r 0 , µ 0,ρ(r) 0

) = W ϕ (µ 1 , ν 1 ) ≤ CW ϕ (µ 0,r * 0 , µ 0,r + 0
) ≤ Cα D (0, r + ). (6.11) In order to show that Tan(µ, 0) is one-dimensional, we define a specific sequence of measures µ 0,r j , which will be used to approximate all the tangent measures of Tan(µ, 0) up to a multiplicative constant. We start with r 0 = 1/4, and define r j by induction, taking r j+1 = ρ(r j ) for j ≥ 0. Note that for all choice of integers 0 ≤ k ≤ l, by (5.6), (6.11), and (6.3) we have

W ϕ (µ 0,r k 0 , µ 0,r l+1 0 ) ≤ k≤j≤l W ϕ (µ 0,r j 0 , µ 0,r j+1 0 ) ≤ C k≤j<l α D (0, (r j ) + ) ≤ 2C k≤j<l ˆ4r j 2r j α D (0, t) dt t .
(6.12)

Recall that r j+1 = ρ(r j ) ≤ λ -1 1 r j (see (6.10)), thus the r j 's decay at a definite rate. Therefore the intervals [2r j , 4r j ] have bounded overlap, and since they are all contained in (0, 4r k ], we obtain

W ϕ (µ 0,r k 0 , µ 0,r l+1 0 ) ≤ C ˆ4r k 0 α D (0, t) dt t . (6.13)
Let σ ∈ Tan(µ, 0) be given. There exist sequences {ρ k } and {a k } such that ρ k ∈ (0, 1/4], lim k→∞ ρ k = 0, a k > 0, and σ k = a k µ 0,ρ k converges weakly to σ.

(6.14) Let j = j(k) denote the largest integer such that r j ≥ ρ k . Thus j ≥ 0 (because r 0 = 1/4), and r j+1 < ρ k ; since r j+1 = ρ(r j ) ∈ [λ -1 2 r j , λ -1 1 r j ] (by (6.10)), we get that

λ -1 2 r j(k) < ρ k ≤ r j(k) . (6.15) Set θ k = ρ k /r j(k) ∈ [λ -1 2 , 1
]; we may replace {ρ k } by a subsequence such that the θ k converge to a limit θ. Consider the dilation D k defined by D k (u) = θ k u, and set D(u) = θu. Notice that by (6.14) and (1.11)

[D k ] ♯ σ k = a k [D k ] ♯ µ 0,ρ k = a k µ 0,ρ k /θ k = a k µ 0,r j(k) . (6.16) Also, the [D k ] ♯ σ k converge weakly to D ♯ σ.
In fact for f continuous and compactly supported,

lim k→∞ ˆf d[D k ] ♯ σ k = lim k→∞ ˆf (θ -1 k x)dσ k (x) = lim k→∞ ˆf (θ -1 x)dσ k (x) = ˆf (θ -1 x)dσ(x) = ˆf d[D] ♯ σ (6.17)
(use the uniform continuity of f and local uniform bounds on the σ k ). By (6.1), σ is a flat measure; then D ♯ σ = σ and (6.17) shows that {a k µ r j(k) } converges weakly to σ. (6.18)

If σ ′ is another nonzero element of Tan(µ, 0), we can find other sequences {j ′ (k)} and {a ′ k }, with lim k→∞ j ′ (k) = ∞ (by the analogue of (6.15)), such that {a ′ k µ 0,r j ′ (k) } converges weakly to σ ′ . (6.19) By Lemma 5.4, then (5.5), and then (6.13) and (6.2),

W ϕ (σ, σ ′ ) ≤ lim inf k→∞ W ϕ (a k µ 0,r j(k) , a ′ k µ 0,r j ′ (k) ) = lim inf k→∞ W ϕ (µ 0,r j(k) , µ 0,r j ′ (k) ) ≤ C lim inf k→∞ ˆ4 max(r j(k) ,r j ′ (k) ) 0 α D (0, t) dt t = 0. (6.20)
Then σ = σ ′ , and this completes our proof of Lemma 6.1.

7 The decomposition of the "good set" in Σ 1

In (1.28)-(1.29) we announced a decomposition of

Σ 0 ∩(Σ 1 ∪Σ 2 ) into pieces S d (0 ≤ d ≤ n),
which satisfy the property that for each x ∈ S d , Tan(µ, x) ⊂ F d . In this section we check that the pieces Σ 1 ∩ S d satisfy the requirements of Theorem 1.2. The remaining sets Σ 2 ∩ S d will be treated in Section 8. We start with d = 0. Set

S 0 = x ∈ Σ ; Tan(µ, x) ⊂ F 0 . (7.1)
We claim (as in the statement of Theorem 1.2) that S 0 is the set of points where µ has an atom, and that every point of S 0 is an isolated point of Σ.

Suppose that µ has an atom at x. Then since µ is doubling, x is an isolated point of Σ (Lemma 2.3 in [START_REF] Azzam | Wasserstein distance and rectifiability of doubling measures: Part I[END_REF]). We can check by hand that Tan(µ, x) is the set F 0 of multiples of the Dirac measure at the origin, and that x ∈ Σ 0 ∩ Σ 1 ∩ Σ 2 (because α D (x, r) = 0 for r small).

Conversely, suppose that Tan(µ, x) ⊂ F 0 , and let us check that x is an isolated point of Σ. Suppose instead that we can find a sequence {x k } in Σ \ {x} that converges to x. Set r k = 2|x -x k |. Since µ is doubling, there is a subsequence of {µ x,r k 0 } which converges weakly to a measure σ. Since σ ∈ Tan(µ, x), σ is a Dirac mass. Let ζ be smooth function such that 1 B\B(0,1/4) ≤ ζ ≤ 1 (B(0,2)\B(0,1/4) . Then ´ζdσ = 0, so lim k→∞ ´ζµ x,r k 0 = 0. On the other hand, by (1.11) and (1.2)

ˆζdµ x,r k 0 = µ(B(x, r k )) -1 ˆζdµ x,r k (7.2) = µ(B(x, r k )) -1 ˆζ(r -1 k (y -x))dµ(y) ≥ µ(B(x, r k )) -1 µ(B(x k , r k /4)) ≥ C -3 δ because ζ(r -1 k (y -x)) = 1 for y ∈ B(x k , r k /4
). This contradiction shows that if x ∈ S 0 , then x is an isolated point in Σ, and then µ has a Dirac mass at x. This gives the desired description of S 0 , the fact that S 0 is at most countable is easy to see.

We may now concentrate on exponents

d ∈ [1, n]. Set S ′ d = x ∈ Σ 0 ∩ Σ 1 ; Tan(µ, x) ⊂ F d = S d ∩ Σ 1 , (7.3)
where the last equality comes from (1.28). Together with S 0 , these sets are disjoint and cover The fact that E = S ′ d satisfies the assumption of the lemma comes from Lemma 6.1.

Σ 1 ∩ Σ 0 (
Proof. If d = n, R n is a d-dimensional
Lipschitz graph that covers E, thus we assume that d < n. We claim that for x ∈ E

x + V x is a tangent plane to Σ at x.

If not there is a sequence {y k } in Σ \ {x}, that tends to x, and such that

dist(y k , x + V x ) ≥ c|y k -x| (7.6)
for some c > 0. Set r k = 2|y k -x|, and replace {y k } with a subsequence for which the {µ x,r k 0 } converges weakly to a measure σ ∈ Tan(µ, x). Let ζ be a smooth compactly supported nonnegative function such that ζ(0) = 0 on V x , but

ζ(u) = 1 for u ∈ B such that dist(u, V x ) ≥ c/2. (7.7)
By assumption, σ is supported on V x and so ´ζdσ = 0. Thus lim k→∞ ´ζµ x,r k 0 = 0. On the other hand, (7.6) says that for y ∈ B(y k , cr k /4), This contradiction proves (7.5).

dist(r -1 k (y -x), V x ) = r -1 k dist(y, x + V x ) ≥ r -1 k [dist(y k , x + V x ) -
For ε > 0 small and x ∈ E, choose an integer j = j(x) ≥ 0 such that dist(y, x + V x ) ≤ ε|y -x| for y ∈ Σ ∩ B(x, 2 -j(x) ). (7.9)

On the Grassmann manifold G(d, n) of vector spaces of dimension d in R n , let us for instance use the distance defined by dist(V, W ) = ||π V -π W ||, where π V and π W denote the orthogonal projections on V and W . With this distance, G(d, n) is compact. Choose a finite family V in G(d, n) such that dist(V, V ) ≤ ε for V ∈ G(d, n). Set E(V, j) = x ∈ E ; j(x) = j and dist(V x , V ) ≤ ε (7.10) for V ∈ V and j ≥ 0. We now cover each E(V, j) with a countable collection of ddimensional Lipschitz graphs. We claim that for each ball B of radius 2 -j-1 , E(V, j) ∩ B is contained in a Lipschitz graph over V . (7.11) Lemma 7.1 follows from this claim because the E(V, j) cover E. To prove the claim, let x, y ∈ E(V, j) ∩ B be given. Observe that |x -y| < 2 -j and y ∈ E ⊂ Σ, so (7.9) guarantees that dist(y, x + V x ) ≤ ε|y -x|. Then

|π V (y) -π V (x)| ≤ |π Vx (y) -π Vx (x)| + ||π V -π Vx |||x -y| ≤ 2ε|x -y|,
which yields (7.11). This completes our proof of Lemma 7.1.

8 The decomposition of the "good set" in Σ 2

Our goal in the section is to apply Theorem 1.5 in [START_REF] Azzam | Wasserstein distance and rectifiability of doubling measures: Part I[END_REF] to the set Σ 0 ∩ Σ 2 , to obtain the desired decomposition. For the reader's convenience we include the necessary background below.

Theorem 8.1 (Theorem 1.5, [START_REF] Azzam | Wasserstein distance and rectifiability of doubling measures: Part I[END_REF]). Let µ be a doubling measure in R n , denote by Σ its support, and set 1. First, Σ 0 (0) is the set of points of Σ where µ has an atom; it is at most countable and each of its point is an isolated point of Σ. We want to apply Theorem 8.1, so we need to show that for each x ∈ Σ 0 ∩ Σ 2 , ˆ1 0 α(x, t) dt t < ∞. (8.1)

Σ 0 = x ∈ Σ ; ˆ1 0 α(x,
V = c V H d V = c V 1 V H d , with c V = H d (V ∩ B)
Let x ∈ Σ 0 ∩ Σ 2 be given. By Theorem 1.3 every tangent measure σ ∈ Tan(µ, x) is flat. To estimate to the distance from µ x,r 0 to σ we proceed as in Section 6 except that we work with the whole group G rather than D. We now follow that argument, without some of the details but we do emphasize the differences.

Without loss of generality, we assume that x = 0. We use the definition of α G and Chebyshev's inequality to associate to each r ∈ (0, 1/4] a radius r + ∈ [2r, 4r] such that α G (0, r + ) ≤ 2 ˆ4r 2r α G (0, t) dt t (8.2)

(see (6.3)). By the definition (1.13), there exists G ∈ G (0, r + ) such that W 1 (µ G 0 , µ 0,r + 0

) ≤ 2α G (0, r + ) (8.3) (see (6.4)). By (1.6) G = λR for some isometry, and since G ∈ G (x), (1.8) guarantee that G(0) = 0 and hence R(0) = 0. That is, R is a linear isometry. We still have that λ ) ≤ 2α G (0, r + ) (8.6)

direction is given by the support of σ; see (7.5). Then S ′′ d is rectifiable, and even satisfies (7.4); see the proof above, and also [START_REF] Azzam | Wasserstein distance and rectifiability of doubling measures: Part I[END_REF]. This completes the proof of Theorem 1.2.

As mentioned in Section 1 we have additional control on the size of µ on Σ 2 and the behavior of µ on S ′′ d . For 1 ≤ d ≤ n and every point x ∈ S ′′ d = Σ 2 ∩ S d , the density of µ exists, that is θ d (x) = lim r→0 r -d µ(B(x, r)) ∈ (0, ∞) (8.18)

(see [START_REF] Azzam | Wasserstein distance and rectifiability of doubling measures: Part I[END_REF]). Morever, we have the further decomposition of S ′′ d into sets S ′′ d ∩ Σ 0 (d, k) where µ and H d are mutually absolutely continuous, as in 4. in Theorem 8.1.

1 7 2

 17 Statement of Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Control of the averages α G 8

  ) = 0 and D(x) = D ∈ D; D(x) = 0 . (1.8)

  .37) Recall that ρ k ≤ t k ≤ 2ρ k by(3.11), and that r -1 k z k is bounded by (3.28); since r -1 k ρ k tends to +∞ by (3.8) we get that for k large, |z k | < Cr k < ρ k ≤ t k , and so B(z k , t k ) ⊂ B(0, 2t k ) ⊂ B(0, 4ρ k ). By (3.31), since 0 ∈ Σ, and by (1.2),

  .39) Return to (3.30), set b k = e k /e ′ k , and observe that by (3.39)

  α k the right-hand side of (3.40) (see (3.22)). Since α k tends to 0 by (3.5), we can choose ρ k and η k so that the constraints (3.8) and (3.6) hold, but the convergence in (3.8) is slow enough so that lim k→∞ β k = 0. (3.41) Recall from (3.34) that C -1 ≤ b k ≤ 1, hence modulo passing to a subsequence (which we relabel) we can guarantee that lim k→∞ b k = b > 0. Letting k tend to infinity in (3.40) we obtain lim k→∞ ˆϕdσ k = ˆϕdσ. (3.42)

  .45) Combining (3.40), (3.41), (3.42), (3.43), (3.44) and (3.45) we obtain that for any 1-Lipschitz function ϕ with compact support, b ˆϕ • H dσ = ˆϕ dσ. (3.46) Since the Radon measure σ is regular, (3.46) also holds for characteristic functions of Borel sets. Hence bH ♯ σ = σ. Recall that λ(H) -1 ∈ [λ 1 , λ 2 ] and H(y) = y by (3.29); thus the conclusion of Lemma 3.1 hold, with c = b -1 .

Lemma 4 . 3 .

 43 .15) Let z k be the orthogonal projection of y k into the line containing x and y. Note that by (4.14) z k ∈ x -2ε|x -y| y-x |x-y| , y + 2ε|x -y| y-x |x-y| . By (4.15), |z k+1 -z k | ≤ 2λ -ℓ 1 |x -y| and by the definition of z k , |z k -y k | ≤ 2ε|x -y|. Therefore every point of [x, y] lies within 2(ε + λ -ℓ 1 )|x -y| of some y k , that is each point of [x, y] is at most 2(ε + λ -ℓ 1 )|x -y| away from Ξ. Choosing ε arbitrarly small and ℓ arbitrarly large, we get that [x, y] ⊂ Ξ. Lemma 4.2 follows. The set Ξ is a vector subspace of R n .

Lemma 4 . 4 .

 44 There is a dimension D ≥ 0 such that σ is Ahlfors-regular of dimension D, which means that C -1 ρ D ≤ σ(B(y, ρ)) ≤ Cρ D for y ∈ Ξ and ρ > 0 (4.18)for some constant C ≥ 1.Proof. Set λ(y) = λ(H(y)), and recall from (4.1) that 1 < λ 1 ≤ λ(y) -1 ≤ λ 2 . Notice that by (4.2), for y ∈ Ξ and r > 0σ(B(y, λ(y) -1 r)) = σ(H(y) -1 (B(y, r))) = H(y) ♯ σ(B(y, r))= c(y)σ(B(y, r)). (4.19) Iterating we obtain for m ≥ 0 σ(B(y, λ(y) -m r)) = c(y) m σ(B(y, r)). (4.20) Applying (4.20) to λ(y) m r instead of r we have that (4.20) also holds for m ≤ 0. Observe that since λ(y) < 1 and σ(B(y, r)) > 0 when y ∈ Ξ, (4.20) yields c(y) ≥ 1. If c(y) = 1, then σ(B(y, λ(y) -m r)) = σ(B(y, r)) for all m ∈ Z, and σ is a Dirac mass. This case was excluded before the statement of the lemma, so c(y) > 1. Now let ρ > 0 be given, and choose m such that λ(y) -m ≤ ρ ≤ λ(y) -m-1 . (4.21) By (4.20) applied to r = 1, c(y) m σ(B(y, 1)) ≤ σ(B(y, ρ)) ≤ c(y) m+1 σ(B(y, 1)), (4.22) hence, letting ℓ = log(σ(B(y, 1))), (4.22) yields m log(c(y)) + ℓ ≤ log(σ(B(y, ρ))) ≤ (m + 1) log(c(y)) + ℓ. (4.23) By (4.21) m log(λ(y) -1 ) ≤ log ρ ≤ (m + 1) log(λ(y) -1 ). (4.24) Hence combining (4.23) and (4.24) we have lim ρ→+∞ log(σ(B(y, ρ))) log ρ = log(c(y)) log(λ(y) -1 ) =: D(y). (4.25) We claim that D(y) does not depend on y. Indeed, if z ∈ Ξ, observe that B(y, ρ) ⊂ B(z, ρ + |z -y|), hence D(y) = lim ρ→+∞ log(σ(B(y, ρ))) log ρ ≤ lim inf ρ→+∞ log(σ(B(z, ρ + |z -y|))) log ρ = lim inf ρ→+∞ log(σ(B(z, ρ + |z -y|))) log(ρ + |z -y|) = D(z); (4.26) the opposite inequality also holds exchanging the roles of y and z. Let D be the common value of the D(y) for y ∈ Ξ. By definition (see (4.25)), λ(y) -D = c(y), and using (4.21) we can rewrite (4.22) as λ(y) -mD σ(B(y, 1)) ≤ σ(B(y, ρ)) ≤ λ(y) -D λ(y) -mD σ(B(y, 1)). (4.27) Thus, by (4.21), λ(y) D ρ D σ(B(y, 1)) ≤ σ(B(y, ρ)) ≤ λ(y) -D ρ D σ(B(y, 1)), (4.28) which yields (4.18) with C = λ D 2 max(σ(B(y, 1)), σ(B(y, 1)) -1 ) (because (4.1) guarantees that λ(y) -1 ≤ λ 2 ).

Lemma 4 . 5 .

 45 Let d be the dimension of the vector space Ξ, and denote by ν = H d Ξ the restriction of H d to Ξ. Let D be as in(4.18). Then d = D, and there exists a constant c 0 > 0 such that σ = c 0 ν.

  .29) and Ξ is a D-dimensional Ahlfors regular set. Hence D = d.By (4.29), σ is absolutely continuous with respect to ν, and the Radon-Nikodym derivative of σ with respect to ν is bounded. Thus there is a bounded function f on Ξ such that σ = f ν.We now show that f is constant. First observe that since d = D = log(c(y)) log(λ(y) -1 ) (by (4.25)), (4.20) guarantees that for y ∈ Ξ, r > 0, and m ∈ Z σ(B(y, λ(y) -m r)) = c(y) m σ(B(y, r)) = λ(y) -md σ(B(y, r)).

  ν(B(y, λ(y) -m r)) = λ(y) -md ν(B(y, r)), we may rewrite (4.30) as σ(B(y, λ(y) -m r)) ν(B(y, λ(y) -m r)) = σ(B(y, r)) ν(B(y, r)) . (4.31) The Lebesgue differentiation theorem says that for ν-almost every y ∈ Ξ, y and every r > 0, by (4.32) and (4.31) for -m f (y) = lim m→∞ σ(B(y, λ(y) m r)) ν(B(y, λ(y) m r)) = σ(B(y, r)) ν(B(y, r)) . (4.33) That is, σ(B(y, r)) = f (y)ν(B(y, r)) for r > 0. (4.34) If z ∈ Ξ is another Lebesgue point of f , since B(z, r) ⊂ B(y, r + |y -z|) we have f (z) = lim r→∞ σ(B(z, r)) ν(B(z, r)) ≤ lim inf r→∞ σ(B(y, r + |y -z|)) ν(B(z, r)) = lim inf r→∞ σ(B(y, r + |y -z|)) ν(B(y, r + |y -z|)) = f (y). (4.35)Similarly f (y) ≤ f (z), and f is constant.

1 ) ; ( 5 . 1 )

 151 If µ and ν are two Radon measures such thatµ(B(0, 1/2)) > 0 and ν(B(0, 1/2)) > 0,(5.2)we define a new distance W ϕ (µ, ν) byW ϕ (µ, ν) = sup ψ∈Lip 1

  c) + c(d -b)| bd (5.12) with a = ´ψϕdµ, b = ´ϕdµ, c = ´ψϕdν, and d = ´ϕdν. Notice that |c| ≤ d because ψ ∈ Lip 1 (B) and ϕ ≥ 0. Also, b ≥ µ(B(0, 1/2)) by (5.1). Hence by (5.9) and (5.10); we have

cr k 4 ] ≥ c/ 4 ,

 44 hence ζ(r -1k (y -x)) = 1 by (7.7), and (1.11) and (1.2) imply thatˆζdµ x,r k 0 = µ(B(x, r k )) -1 ˆζdµ x,r k (7.8) = µ(B(x, r k )) -1 ˆζ(r -1 k (y -x))dµ(y) ≥ µ(B(x, r k )) -1 µ(B(x k , cr k /4)) ≥ C -1 .

  , r) = min d=0,1,...,n α d (x, r), and α d (x, r) = inf W 1 (µ x,r 0 , ν V ) ; V ∈ A ′ (d, n) .Here A ′ (d, n) denotes the set of n dimensional affine spaces which intersect B(0, 1/2) and ν

- 1 .

 1 Then there are disjoint Borel setsΣ 0 (d) ⊂ Σ, 0 ≤ d ≤ n, such that Σ 0 = n d=0 Σ 0 (d),with the following properties.

2 .

 2 For 1 ≤ d ≤ n and x ∈ Σ 0 (d), the limitθ d (x) := lim r→0 r -d µ(B(x, r))exists, and 0 < θ d (x) < ∞.

3 . 4 .

 34 For 1 ≤ d ≤ n and x ∈ Σ 0 (d), Σ has a tangent d-plane at x, W , and set W * = W -x. Then Tan(x, µ) = {cH d W * ; c > 0}. In addition, the measures µ x,r 0 converge weakly to H d W * . Further decompose Σ 0 (d), 1 ≤ d ≤ n, into the setsΣ 0 (d, k) = x ∈ Σ 0 (d) ; 2 k ≤ θ d (x) < 2 k+1 , k ∈ Z; then each Σ 0 (d, k) is a rectifiable set of dimension d, with H d (Σ 0 (d, k) ∩ B(0, R)) < ∞ for every R > 0,µ and H d are mutually absolutely continuous on Σ 0 (d, k), and µ = θ d H d there.

- 1 =

 1 λ(G) -1 ∈ [λ 1 r + , λ 2 r + ], and if we setr * = λ -1 ∈ [λ 1 r + , λ 2 r + ]; (8.4)as in (6.5), we have thatG(B(0, r * )) = B (8.5)and µ G 0 is the image of µ 0,r * 0 by a linear isometry. That is, µ G 0 = R ♯ µ 0,r * 0 and (8.3) only yields W 1 (R ♯ µ 0,r * 0 , µ 0,r + 0

  that depends on C δ and λ 2 . To estimate a k , consider a test function f such that 1 B ≤ f ≤ 1 2B . By definition of σ, ´f dσ = lim k→∞ ´f dσ k . By (1.2), (3.3) and the definition above (1.11), we have

	by (1.12). Then (1.2), (3.31), (3.32), (3.33), and the fact that λ 1 > 1
	yield
	C -1 e ′ k ≤ e k ≤ e ′ k	(3.34)
	for some constant C

  by (1.29)), hence also µ-almost all of Σ 1 (by (1.26)). By Lemma 6.1, the only part of Theorem 1.2 concerning Σ 1 that remains to be proven is the fact that S ′ d is rectifiable for 1 ≤ d ≤ n, and more precisely Suppose that µ is a doubling measure, Σ is its support, d ∈ {1, ..., n}, and E ⊂ Σ is such that for all x ∈ E, there is a vector space V x of dimension d such that Tan(x, µ) = cH d | Vx ; c > 0 . Then E can be covered by a countable family of d-dimensional Lipschitz graphs.

	S ′ d can be covered by a countable family	(7.4)
	of d-dimensional Lipschitz graphs.
	(This is slightly more precise because we don not need to add a H d -negligible set.) This
	follows from the following lemma, which is essentially known, but which we prove for the
	reader's convenience.
	Lemma 7.1.

Acknowledgements

The authors are grateful to Alessio Figalli and Xavier Tolsa for helpful discussions. Parts of this work were done while the first author was visiting IPAM and while the second author was visiting the University of Washington.

The first author was partially supported by NSF RTG grant 0838212. The second author acknowledges the generous support of the Institut Universitaire de France, and of the ANR (programme blanc GEOMETRYA, ANR-12-BS01-0014). The third author was partially supported by an NSF grants DMS-0856687 and DMS-1361823, a grant from the Simons Foundation (# 228118) and the Robert R. & Elaine F. Phelps Professorship in Mathematics.

instead of (6.6). We still multiply the radii by r/r * , set

as in (6.10), and deduce from (8.6) that W ϕ (R ♯ µ 0,r 0 , µ 0,ρ(r) 0

) ≤ Cα G (0, r + ), (8.8) using the same proof which involves Lemma 5.1 and Lemma 5.2 (the extra rotation does not affect the argument). Inequality (8.8) is the analogue of (6.11). Let us write this slightly differently. Set R r = R -1 ; then by (8.8)

since the W ϕ -distance is invariant under isometry.

Given r 0 ≤ 1/4, we can construct a decreasing sequence {r j } as we did before, defined by r j+1 = ρ(r j ). Let us keep track of the rotations: set S 0 = I and S j+1 = S j R r j . For k ≥ 0 we want to estimate the numbers

Let us check by induction that

When k = 0, this is (8.9) for r 0 . If k ≥ 1 and (8.11) holds for k -1, the triangle inequality (5.6) yields

by definition of S k+1 , the invariance of W ϕ under linear isometries, and (8.9). This proves (8.11). Then (8.2) and the same argument as in (6.12)-(6.13) yield

The final integral is finite because 0 ∈ Σ 2 (see the definition (1.17)).

The measures µ 0,r k 0 are suitably normalized, so there is a subsequence which converges weakly to some measure σ (again see Lemma 2.1 in [START_REF] Azzam | Wasserstein distance and rectifiability of doubling measures: Part I[END_REF] for a little more detail). There is also a further subsequence for which the S k converge to an isometry S, and then the S k ♯ µ 0,r k 0 converge to S ♯ σ (proceed as for (6.17)). By Lemma 5.4, (8.10), and (8.13),

.14)

We now use Lemma 5.3 to translate estimate (8.14) into an upper bound for the ´1 0 α(x, r) dr r . For t ∈ [1/4, 1/2], the measure µ t that is defined by (5.27) with µ replaced by µ 0,r 0 0 is just µ 0,tr 0 0 . Since σ is a flat measure, so is S ♯ σ. Hence the measure ν t built from ν = S ♯ σ as in (5.27) is also a flat measure supported on a d-plane V passing through the origin. We use ν t to estimate α(x, r). By (5.28), the fact that µ is doubling, and (8.14), we have

Note that (8.15) holds for r 0 ≤ 1/4.

We are now ready to prove that for x ∈ Σ 0 ∩ Σ 2 , ´1 0 α(x, r) dr r < ∞. Recall we are assuming x = 0. By (8.15) and the definition (1.17) of Σ 2 we have

A brutal estimate shows that α(x, r) ≤ 2 for r > 0, hence ´1 1/8 α(x, s) ds s < ∞. Hence the hypothesis of Theorem 8.1 hold. We obtain a decomposition of the set Σ 0 ∩ Σ 2 into subsets S ′′ d , 0 ≤ d ≤ n, which satisfy all the requirements for Theorem 1.2. In fact we get some additional information which we record here. First, for every point x ∈ S ′′ d , d ≥ 1, Tan(µ, x) is the vector space of dimension 1 spanned by some flat measure of dimension d (see (1.27)) which implies S ′′ d = x ∈ Σ 0 ∩ Σ 2 ; Tan(µ, x) ⊂ F d = Σ 2 ∩ S d , (8.17)

where S d is as in (1.28). Moreover once we know that Tan(µ, x) is the space of dimension 1 spanned by some flat measure σ ∈ F d , we have that Σ has a tangent d-plane at x whose


P R Arantes Gilz

C Louembet

Predictive control algorithm for spacecraft rendezvous hovering phases

Keywords: predictive control, rendezvous, periodic relative motion, constrained trajectories, non-negative polynomials, alternating projections algorithm

The present article proposes a predictive control law for orbital rendezvous hovering phases. The control consists of impulsive thrust corrections that maintain the relative motion between spacecraft bounded. The proposed predictive algorithm accounts for thrusters saturation and a propellant budget on top of the state constraints. It is shown that computing an admissible control is equivalent to determine a point in the intersection between two closed convex sets. To resolve this problem, an alternating projections algorithm is used. The efficiency of the proposed predictive control algorithm is bench tested on a nonlinear simulator.

I. INTRODUCTION

In the future missions of on-orbit inspection and on-orbit servicing, the ability to maintain a relative periodic motion is a crucial issue for any rendezvous spacecraft. The purpose of this research is to develop a predictive control law that enables the control of periodic relative orbits that are included in a given linear region for sake of safety.

Constrained relative trajectories have been studied under the term of "hovering". Hovering has been defined in [START_REF] Irvin | Fuel-optimal maneuvers for constrained relative satellite orbits[END_REF] as follow

To develop a control strategy to place a deputy satellite inside a specific lobe defined in the chief body-fixed frame and keep it there in a fuel-optimal manner.

In most of works on the subject, as [START_REF] Irvin | Fuel-optimal maneuvers for constrained relative satellite orbits[END_REF]- [START_REF] Williams | Optimal impulsive maneuvering within a confined hover region[END_REF], the hovering strategy is based on two different strategies. First, the "teardrop" strategy consists in computing an impulsive thrust such that the initial and final relative positions on the bounds of the tolerance subset are equal for a given time of flight [START_REF] Irvin | Fuel-optimal maneuvers for constrained relative satellite orbits[END_REF], [START_REF] Irvin | An investigation of teardrop relative orbits for circular and elliptical chief satellites[END_REF]. Second, the "pogo" strategy, where the impulsive control is computed each time the vehicle hit the bounds of the lobe to maintain it as longtime as possible while minimizing the fuel consumption [START_REF] Irvin | Fuel-optimal maneuvers for constrained relative satellite orbits[END_REF], [START_REF] Luo | Optimal multi-objective linearized impulsive rendezvous[END_REF], [START_REF] Williams | Optimal impulsive maneuvering within a confined hover region[END_REF].

However, some contributions have recently taken advantage of the naturally periodic relative orbits to design hovering control laws [START_REF] Deaconu | Constrained periodic spacecraft relative motion using non-negative polynomials[END_REF]- [START_REF] Deaconu | Designing continuously constrained spacecraft relative trajectories for proximity operations[END_REF]. In [START_REF] Deaconu | Constrained periodic spacecraft relative motion using non-negative polynomials[END_REF], a methodology to compute the coordinates of periodic relative orbits included in a polytopic region is exposed. This method has been extended to the non periodic relative trajectories in [START_REF] Deaconu | Designing continuously constrained spacecraft relative trajectories for proximity operations[END_REF].

A control law proposed in [START_REF] Deaconu | On the trajectory design, guidance and control for spacecraft rendezvous and proximity operations[END_REF] consists of two impulsive maneuvers that stabilize the spacecraft with respect to a given constrained periodic relative orbit previously computed. Even though this control law has been proven stable for any eccentric target orbit and requires very few computational effort, its optimality in terms of fuel consumption is not guaranteed.

The present article proposes a different strategy to hover inside a given polytopic subset of the relative state space. The purpose is to steer the vehicle to the set of the periodic relative orbits enclosed in a hovering zone.

To achieve such a goal, a predictive methodology is applied to provide the impulsive control. At each call, the controller computes, if it exists, the coordinates of a relative periodic orbit included in the hovering zone that go through the current position. The impulsive control is then deduced from these coordinates. Moreover, during the control computation, spacecraft abilities and propellant resource are also taken into account.

Although the present work extends the results from the conference article [START_REF] Deaconu | A two-impulse method for stabilizing the spacecraft relative motion with respect to a periodic trajectory[END_REF], it differs in several points. In the present paper, the chaser is steered to a periodic and admissible orbit not defined a priori, contrary to what is done in [START_REF] Deaconu | A two-impulse method for stabilizing the spacecraft relative motion with respect to a periodic trajectory[END_REF]. Consequently, the control consists of one impulse instead of two for the techniques from [START_REF] Deaconu | A two-impulse method for stabilizing the spacecraft relative motion with respect to a periodic trajectory[END_REF]. Finally, the proposed methodology permits to consider in the problem saturation and propellant budget constraints.

As it is exposed thereafter, the set of the orbits coordinates admissible with respect to the above-mentioned constraints is described as the intersection between a semi-algebraic set and a hyperplane. Obtaining the control mainly consists in computing one point of this intersection set. This is achieved using an alternating projections algorithm. Note that this algorithm by means of its properties make the use of a complete optimization procedure unnecessary as explained in this paper.

The paper is structured as follows. Section II gives a brief description of the relative motion and recall the trajectories parametrization used in the sequel. Section III describes, in terms of the relative trajectories coordinates, the constraints on the relative orbits and the controls. Section IV presents the alternating projection algorithm that permits to compute the coordinates of an admissible relative trajectory and the control. Finally, section VI proposes numerical examples.

II. RELATIVE MOTION MODELING

The Fig. 1 presents the frames used to model the relative motion between the leader S l and the follower S f spacecraft. suming Keplerian orbits, the relative dynamics can be given by the linearized and simplified Tschauner-Hempel equations [START_REF] Tschauner | Elliptic orbit rendezvous[END_REF] (time is replaced by true anomaly as independent variable):

x′′ = 2z ′ + ũx ỹ′′ = -ỹ ′ + ũy z′′ = 3 1+e cos ν z -2x ′ + ũz (1)
where e is the leader orbit eccentricity. Defining the vector X(ν) = [x(ν) ỹ(ν) z(ν) x′ (ν) ỹ′ (ν) z′ (ν)] T , a state-space representation for the non-autonomous system (1) is:

X′ (ν) = Ã(ν) X(ν) + B ũ (2)
Note that the state expressed in time and expressed in true anomaly are linked by means of the following similar transformation:

X(ν) = (1 + e cos ν)I 3 0 3 -e sin νI 3 √ a 3 (1-e 2 ) 3 √ µ(1+e cos ν) 2 I 3 X(t) (3)
Using the transition matrix, Φ(ν, ν 0 ), proposed by Yamanka and Ankersen in [START_REF] Yamanaka | New state transition matrix for relative motion on an arbitrary elliptical orbit[END_REF] to propagate a initial state X(ν 0 ) and considering an unforced motion:

X(v) = Φ(ν, ν 0 ) X(ν 0 ), ν ≥ ν 0 (4)
Expanding eq. ( 4), the following relative trajectories description can be deduced:

x(ν) = (2 + e cν )(d1 sν -d2 cν ) + d3 + 3 d0 J ν ν 0 (1 + e cν ) 2 ỹ(ν) = d4 cν + d5 sν z(ν) = (1 + e cν )(d2 sν + d1 cν -3 e d0 J ν ν 0 sν ) + 2 d0 (5)
where c ν = cos ν, s ν = sin ν, J ν ν0 is given by:

J ν ν0 = ν ν0 dτ (1 + e cos τ ) 2 (6)
and

D(ν) = [d 0 (ν) d 1 (ν) d 2 (ν) d 3 (ν) d 4 (ν) d 5 (ν)
] T , the parameter vector, is expressed as a linear function of X(ν),

D(ν) = C(ν) X(ν) (7)
Please refer to [START_REF] Deaconu | Designing continuously constrained spacecraft relative trajectories for proximity operations[END_REF] for a more detailed equation ( 7). In this section, the admissible trajectories are described. The admissible orbits must have two main properties: the periodicity and the inclusion in a polytopic set. These two properties are described thereafter as constraints on the parameter vector D.

A. Periodicity

The relative trajectories between spacecraft are not generally periodic. The periodicity constraint guarantees a relative motion that does not require any correction in absence of disturbances. Imposing X(ν + 2π) = X(ν), ∀ν ≥ ν 0 one can deduce that this constraint can be expressed as:

d 0 (ν 0 ) = 0 (8)
Analyzing the equations in ( 5), condition (8) becomes evident, since the only non-periodic and divergent term J ν ν0 is always multiplied by d 0 .

B. Polytope constraints using non-negative polynomials

During the rendezvous hovering phases, the follower is required to remain in the interior of a certain limited region of the space. This region is going to be modeled as a polytope (without any generality loss, a rectangular parallelepiped will be used):

x min ≤ x(t) ≤ x max y min ≤ y(t) ≤ y max z min ≤ z(t) ≤ z max , ∀t ≥ t 0 (9)
Defining X pos (t) = [x(t) y(t) z(t)] T , this constraint can be rewritten using the matrix H = [I 3×3 -I 3×3 ] T and the vector K = [x max y max z max x min y min z min ] T :

H X pos (t) ≤ K, ∀t ≥ t 0 (10)
Replacing the time t by the true anomaly of the leader ν as independent variable, we obtain:

H Xpos (ν) ≤ K(ν), ∀ν ≥ ν 0 (11)
Assuming that the periodicity constraint is satisfied and introducing the following variable change:

w = tan ν 2 , cos ν = 1 -w 2 1 + w 2 , sin ν = 2w 1 + w 2 (12)
the trigonometric terms in ( 5) are eliminated and we obtain the following rational expressions:

x(w) = 1 (1+w 2 ) 2 4 i=0 p xi w i = 1 (1+w 2 ) 2 P x (w) ỹ(w) = 1 (1+w 2 ) 2 i=0 p yi w i = 1 (1+w 2 ) P y (w) z(w) = 1 (1+w 2 ) 2 4 i=0 p zi w i = 1 (1+w 2 ) 2 P z (w) (13)
the coefficients p xi , p yi and p zi depending linearly on D:

p x0 = (-2 -e)d 2 + d 3 p x1 = (4 + 2e)d 1 p x2 = (2 + e)d 2 + 2d 3 p x3 = (4 -2e)d 1 p x4 = (2 -e)d 2 + d 3 p y0 = d 4 p y1 = 2d 5 p y2 = -d 4 p z0 = (e + 1)d 1 p z1 = (2e + 2)d 2 p z2 = -2ed 1 p z3 = (2 -2e)d 2 p z4 = (e -1)d 1 ( 14
)
Using the variable change [START_REF] Deaconu | A two-impulse method for stabilizing the spacecraft relative motion with respect to a periodic trajectory[END_REF] and the rational expressions ( 13) in [START_REF] Deaconu | Designing continuously constrained spacecraft relative trajectories for proximity operations[END_REF], we can obtain the following polynomial inequalities: P x (w) ≤ ((1e)w 4 + 2w , ∀w ≥ w 0 [START_REF] Nesterov | Squared functional systems and optimization problems[END_REF] that can be reformulated as:

Γ i (w, D) ≥ 0, ∀w ≥ w 0 , i = 1..6 (16)
The non-negativity conditions in [START_REF] Boyd | Alternating projections[END_REF] result in an infinite number of constraints to be respected. Deaconu presents in [START_REF] Deaconu | Constrained periodic spacecraft relative motion using non-negative polynomials[END_REF] a way to transform this infinite number of constraints in a finite number of convex constraints using the following theorem presented by Nesterov in [START_REF] Nesterov | Squared functional systems and optimization problems[END_REF]: Theorem 1: Let P (w) be a polynomial of degree 2n on R and p ∈ R 2n+1 the vector containing its coefficients. Then:

P (w) ≥ 0, ∀w ∈ R ⇔ ∃Y ∈ S n+1 + |p = Λ * (Y ) p l = Λ * l (Y ) = trace(Y H n,l ), l = 1..2n + 1 (17)
where S n+1 + is the cone of the symmetric semi-definite positive matrices of size n + 1 and H n,l are the Henkel matrices:

H n,l (j, k) = 1, if j + k = l + 1 0, otherwise (18)
This is equivalent to affirm that the cone of non-negative univariate polynomials is the image of the cone of positive semi-definite matrices by the linear operator Λ * (•).

Let γ i be the vectors containing the coefficients of the polynomials Γ i (w). Using [START_REF] Cheney | Proximity maps for convex sets[END_REF], we can transform conditions in [START_REF] Boyd | Alternating projections[END_REF] into:

H i Xpos (ν, D) ≤ Ki (ν) Γ i (w, D) ≥ 0, ∀w ≥ w 0 ∃Y i ∈ S + such that γ i = Λ * (Y i ) , i = 1..6 (19)
In the end of this section, constraints on the control are given: the thrusters saturation and the propellant budget. Considering saturation constraints is one of the main advantages of the predictive control law. Moreover, the propellant is a rare resource for spacecraft systems since they can not be refueled. Thus, propellant budget constraints is also proposed to take into account the scarceness of the propellant.

C. Thrusters saturation

Let ∆V = [∆V x ∆V y ∆V z ] T be the vector representing the impulsive velocity correction executed by the follower. Considering that the spacecraft have a pair of equivalent propellers symmetrically and oppositely disposed on each axis and assuming that the saturation limit for each propeller is ∆V max > 0, the saturation constraint can be written as:

|∆V s | ≤ ∆V max , s = x, y, z (20)
These inequalities can be transformed in a matrix positive semi-definite constraint:

|∆V s | ≤ ∆V max ⇔ ∆Vmax ∆Vs ∆Vs ∆Vmax Mi 0, s = x, y, z (21)

D. Propellant budget

An maximum budget of propellant σ is available for each impulsion executed. This constraint is expressed as:

||∆V || 1 = |∆V x | + |∆V y | + |∆V z | ≤ σ (22)
To check condition (22), we must verify if:

∃ Z x , Z y , Z z ≥ 0 W + x , W + y , W + z ≥ 0 W - x , W - y , W - z ≥ 0 W 0 ≥ 0 s.t. For s = x, y, z W - s = Z s -∆V s W + s = Z s + ∆V s σ -s Z s = W 0 (23)

IV. COMPUTING THE CONTROL LAW

Accordingly to the model predictive control scheme, the control is computed periodically in order to bring the system on a periodic trajectory inside a given polytope respect to actuator saturation and budget constraints. Using the previous equations, links between the variables of the problem and its restrictions can be stated. The Fig. 2 presents these links, each ellipse representing a constraint. Given a set of constraints {H, K, ∆V max , σ} and an initial relative state X 0 , we want to determine an impulse ∆V respecting both saturation and budget constraints such that the state after the impulse, X + = X 0 + B∆ Ṽ ,equivalent to the state D + (by means of ( 7)) whose first component d 0 is zero to copt the periodicity and that generates (by means of ( 14), ( 15) and ( 16)) six non-negative polynomials to ensure the polytopic inclusion of the orbits

A. Closed convex sets intersection problem

The constraints can be summarized as follows:

                                               Periodicity: d 0 = 0
Polytope constraints:

∃Y i 0 s.t. γ i = Λ * (Y i ), i = 1..6
Thrusters saturation:

∃ Ms 0 s.t. ∆Vmax ∆Vs ∆Vs ∆Vmax = Ms, s = x, y, z
Propellant budget:

∃ Zs ≥ 0 W + s ≥ 0 W - s ≥ 0 W 0 ≥ 0 s.t. W - s = Zs -∆Vs W + s = Zs + ∆Vs σ -s Zs = W 0 , s = x, y, z (24)
We remark that these constraints are convex, since they are mathematically represented by matrices positive semidefiniteness (on the left) and equalities (on the right).

Let Q be the symmetric block diagonal matrix given by:

Q = diag(M x , M y , M z , Z x , Z y , Z z , W + x , W + y , W + z , W - x , W - y , W - z , W 0 , Y 1 , . . . , Y 6 ) (25)
This matrix is sparse by construction. Table I shows that taking symmetry into account we only need 49 variables to completely describe it. We can use Q to reformulate the positive semi-definiteness conditions and the equations linking the variables in (24) as follows:

tr(A i Q) = b i , A i ∈ S 32 , b i ∈ R, i = 1..m Q 0 (26)
where S 32 is the set of 32 × 32 symmetric matrices and m is the number of equations.

Equations tr(A i Q) = b i describe an affine subset of S 32 and Q 0 represents the symmetric semi-definite matrices cone S 32 + . Therefore, solving the problem is equivalent to find an element Q belonging to the intersection of these sets.

B. Alternating projections algorithm

For two closed and convex sets, we can apply the alternating projections algorithm proposed by Boyd and Dattorro in [START_REF] Boyd | Alternating projections[END_REF] to determine the existence of an intersection: Let C and D be two closed convex subsets of E, let || • || E be a norm on E and let P C and P D denote the projections on C and D, respectively. The algorithm consists in:

Algorithm 1 Alternating Projections Choose an arbitrary Q 0 and compute Q C 1 = P C (Q 0 ) for k = 0, 1, 2, . . . do Q D k = P D (Q C k ) Q C k+1 = P C (Q D k ) end for
It has been proven that [START_REF] Cheney | Proximity maps for convex sets[END_REF]:

• if C ∩ D = ∅, ∃ Q * ∈ C ∩ D s.t. Q C k k→∞ -→ Q * and Q D k k→∞ -→ Q * ; • if C ∩ D = ∅, ∃ Q C * ∈ C, Q D * ∈ D s.t. Q C k k→∞ -→ Q C * and Q D k k→∞ -→ Q D * and ||Q C * -Q D * || E = dist(C, D),
i.e. a pair of points in C and D that have minimum distance.

Thereafter, let C be the affine set described by the equations from system (26), let D be the symmetric semi-definite matrices cone S 32 + and let ||•|| F be the Frobenius norm. Since these sets are closed and convex, the alternating projections algorithm can be implemented to determine an element belonging to their intersection. In this case, projection P C is given by:

P C (Q) = Q - m i=1 a i A i ( 27
)
where a i are solution of the following equation:

Ga = g ( 28
)
where

G(j, k) = trace(A j A k ) g = [trace(A 1 Q) -b 1 . . . trace(A m Q) -b m ] T (29)
The projection P D is given by:

P D (Q) = U diag(max {0, λ 1 } , . . . , max {0, λ 32 }) U T (30)
by considering the spectral decomposition is written as

Q = U diag(λ 1 , . . . , λ 32 ) U T (31)

C. Initialization and stopping criteria

As any iterative algorithm, the projection algorithm's behavior depends on two items: the initial guess and the stopping criteria. The initial guess to provide to the alternating projection algorithm is the matrix Q 0 . In fact, the alternating projection algorithm provides to the initial guess the nearest element in the intersection of the subsets (if exists) [START_REF] Dykstra | An algorithm for restricted least squares regression[END_REF]. Thus, the choice of Q 0 can have a great impact the solution Q * as can be seen of 3. In this paper, two strategies are tested: the initial guess Q 0 is either the null matrix or the solution Q * of the previous call of the projection algorithm. As it is exposed in section VI, different strategies for choosing the initial iterate will lead to different controlled trajectories.

Q 0 C D Q ′ 0 Fig. 3: Alternating projection algorithm
The stopping criteria is based on the distance in the Frobenius norm between

Q D k and Q C k : δ k = log Q D k -Q C k F (32)
The algorithm stops when

δ k ≤ ǫ Q or δ k+1 -δ k ≤ ǫ δ with δ k+1 > ǫ Q (33)
where ǫ Q > 0 and ǫ δ > 0. The first condition means the iterate Q k is close enough to the intersection and the second condition reveals that the algorithm reaches a minimal and non zero distance between the convex subsets so that the intersection is an empty set.

V. STABILITY BY MEANS OF INVARIANCE

The stability around a periodic trajectory is understood in the sense given in [START_REF] Khalil | Nonlinear Systems[END_REF], which requires the characterization of the evolution of the trajectories that start arbitrarily close to the periodic one. According to [20, theorem 3.5], the stability of periodic system is characterized by the eigenvalues of its monodromy matrix M = Φ(ν + 2π, ν). In the case of the linearized relative motion, M has at least two eigenvalues µ i such that |µ i | > 1. This leads to the conclusion that the autonomous system X′ (ν) = Ã(ν) X(ν) is naturally unstable. Thus, without control, the trajectories, described by [START_REF] Luo | Optimal multi-objective linearized impulsive rendezvous[END_REF], that started arbitrary close to any periodic trajectory will diverge from it.

The proposed predictive control law has for goal to steer the system in the vicinity of the set of admissible periodic orbits. Thus, the stability will be defined in the sequel as the invariance of given subsets described thereafter. Now, two sets are given in order to define the stability of the system controlled by means of the proposed predictive algorithm:

• S 0 {H, K}, the set of vectors D(ν) ∈ R 6 that describes periodic trajectories satisfying the polytope constraints;

• S 1 {H, K, ∆V max , σ}, the set of vectors D(ν) ∈ R 6 for which it is possible to find one impulse ∆V , respecting the saturation constraints and the propellant budget, such that D + (ν), the state right after the impulse, belongs to S 0 {H, K}.

The set S 0 {H, K} naturally invariant because of the periodicity respected by each of its elements.

The set S 1 {H, K, ∆V max , σ} is also a class of equivalence under the developed control law, since it is possible to compute an impulsive correction via the alternating projections algorithm that generates a new state belonging to S 0 {H, K}. Consequently, the set S 0 is invariant and attractive under the developed control law with S 1 as convergence set. If the current state does not belong to S 1 , an ad hoc control law must be used to bring the system in S 1 . For instance, an optimal impulsive plan computed using techniques from [START_REF] Deaconu | Designing continuously constrained spacecraft relative trajectories for proximity operations[END_REF] can bring the chaser satellite to the vicinity of the hovering box from any relative state.

The Figure 4 represents the sets S 0 and S 1 and the stability and instability zones in R 6 .

Unstable zone

Attractive zone to S 1 S 1 S 0 Fig. 4 This mission has been implemented on a simulator based on Gauss nonlinear equations for the relative motion, taking into account the atmospheric drag and Earth's oblateness effects (see [START_REF] Vallado | Fundamentals of astrodynamics and applications[END_REF] for details).

In order to assess the impact of the initialization techniques described in IV-C on the resultant trajectory and the consumption, a warm start initialization, taking the last iteration solution as initial guess, and a cold start initialization, using a constant Q 0 matrix for all the iterations, were performed. Simulations with different eccentricities and time-intervals between control calls are conducted to test the algorithm with both initialization procedures. The aim of these simulations is first to verify the designed control law efficiency in respecting the different constraints while rejecting the unmodeled dynamics, and second to highlight the importance of the algorithm initialization procedure. The results obtained for these simulations are presented in Figure 5. We can remark that, on one hand, the smaller eccentricity, the lower consumption is. On the other hand, the consumption increases with the period between controls. Another conclusion can be drawn from Figure 5: the control algorithm using a cold start technique generates greater consumptions compared to the one using a warm start technique. This can be explained by the fact that cold start leads to computation of parameters of admissible orbit that evolve much more than the one computed from warm start procedure.

On Figure 6, controlled trajectories obtained by means of the nonlinear simulator previously described are exposed. The control algorithm is called every 20 second for this specific simulation. First, the simulated trajectories remains admissible during the 5-orbit period of hovering. However, it can be observed that both initialization procedures lead to different geometries respecting the tolerance box constraints despite the orbital perturbation. The total ∆V consumption over the 5 orbital periods is of 4.4 mm/s and 18 mm/s for the warm start and cold start control respectively. Note that, the ∆V budget and saturation constraints are also satisfied. These facts highlight that the admissible orbit computed using a warm start initialization is closer to current position than the one computed from a cold start. In fact, after every call of the control algorithm, the chaser will evolve close to the computed orbit. With the warm start initialization, the next computed orbit will remains in the vicinity of the previously computed orbit. In the contrary, the cold start procedure leads to the computation of a completely different admissible orbit.

VII. CONCLUSION

The present paper develops a predictive control law for the orbital rendezvous. This control law aims to make the chaser satellite hover in a given tolerance set in order to maintain a way point situation. As the relative motion has a periodic dynamics, the predictive control algorithm take advantage of the existence of periodic orbits for this system to develop a law that maintain the hovering while comsuming propellant essentially to reject perturbation, contrary to pogo or tear drop strategies available in the literature. This control law consists in steering periodically the chaser spacecraft into the subset of the admissible orbits. Moreover, budget and saturation condition are also accounted in the control computation. One of the main features of the proposed control law is that it is based on the alternating projection algorithm. This iterative algorithm provides to any initial guess the nearest point in the intersection subset of convex sets when it is not empty. We studied the influence of the initial guess choice on the controlled trajectories and the consumption through nonlinear simulations. Further studies should also focus on the numerical issues like convergence time as function of the stopping criteria tolerance. Another problem of interest is the characterization of the invariant set S 1 , that could speed up the decision of whether or not to apply the projection algorithm, once it should be possible to determine the algorithm convergence by verifying the belonging of the current state to S 1 .

The ECI 1 Fig. 1 :

 11 Fig. 1: Inertial and relative frames

Fig. 2 :

 2 Fig. 2: Variables and its constraints

  -20 20 -20 20] T Hovering duration Time between impulsions [s] 5 orbital periods {20, 50, 100, 200, 500, 1000, 2000}

Fig. 5 :

 5 Fig. 5: Mission simulation results

Fig. 6 :

 6 Fig. 6: Nonlinear simulation for e = 0.1 over 5 orbital periods

  2 + 1 + e)x max P y (w) ≤ y max P z (w) ≤ ((1e)w 4 + 2w 2 + 1 + e)z max P x (w) ≥ ((1e)w 4 + 2w 2 + 1 + e)x min P y (w) ≥ y min P

z (w) ≥ ((1e)w 4 + 2w 2 + 1 + e)z min

TABLE I :

 I Q matrix

		#	Size	Total Size	#Vars	Total #Vars
	Ms	3	2 × 2	6 × 6	3	9
	Zs, W + s , W -s , W 0	10	1 × 1	10 × 10	1	10
	Y {1,3,4,6}	4	3 × 3	12 × 12	6	24
	Y {2,5}	2	2 × 2	4 × 4	3	6
	Total			32 × 32		49

  : Stability zone VI. NUMERICAL EXAMPLES To illustrate our methodology, the example of a given mission is proposed in table II.

TABLE II :

 II Test case mission

	Semi major axis [km]	Eccentricities
	20000	{0.1, 0.2, 0.3, 0.4, 0.5, 0.6}
	Budget [m/s]	Saturations ∆Vmax [m/s]
	0.3	2
	Initial state	Tol. box [xmin, xmax, ymin, ymax, zmin, zmax]