Rodolphe Garbit
email: rodolphe.garbit@univ-angers.fr

A CENTRAL LIMIT THEOREM FOR TWO-DIMENSIONAL RANDOM WALKS IN A CONE

Keywords: . 2000 Mathematics Subject Classification. 60F17, 60G50, 60J05, 60J65 conditioned random walks, Brownian motion, Brownian meander, cone, functional limit theorem, regularly varying sequences

published or not. The documents may come

A central limit theorem for two-dimensional random walks in a cone

1. Introduction 1.1. Main result. The aim of this paper is to underscore a natural necessary and sufficient condition for the weak convergence of a two-dimensional random walk conditioned to stay in a cone to the corresponding Brownian meander. The condition only involves the asymptotic behavior of the tail distribution of the first exit time from the cone.

Let (ξ n ) n≥1 be a sequence of independent and identically distributed random vectors of R d , d ≥ 1, defined on a probability space (Ω, F, P). We assume that the distribution of ξ 1 satisfies E(ξ 1 ) = 0 and Cov(ξ 1 ) = σ 2 I d , where σ 2 > 0 and I d is the d × d identity matrix.

We form the random walk S = (S n ) n≥1 by setting S n = ξ 1 + • • • + ξ n , and for each n ≥ 1, we define a normed and linearly interpolated version of S by

S n (t) = S [nt] σ √ n + (nt -[nt]) ξ [nt]+1 σ √ n , t ≥ 0,
where [a] denotes the integer part of a.

The weak convergence of the process S n = (S n (t), t ≥ 0) as n → ∞ to a standard Brownian motion is Donsker's theorem (see for example Theorem 10.1 of [START_REF] Billingsley | Convergence of Probability Measures[END_REF]).

We consider a linear cone C ⊂ R d (i.e. λC = C for every λ > 0) with the following properties:

(1) C is convex, (2) its interior C o is non-empty,

P(ξ 1 ∈ C \ {0}) > 0. (3)
Such a cone is said to be adapted to the random walk. Note that the convexity of C ensures that its boundary ∂C is negligible with respect to Lebesgue measure (see for example [START_REF] Lang | A note on the measurability of convex sets[END_REF]). The third condition ensures that the first step of the random walk is in C with positive probability. Since a convex cone is a semi-group, the event {ξ 1 , . . . , ξ n ∈ C} is a subset of {S 1 , . . . , S n ∈ C}, so the latter has also a positive probability. For this purpose, one could simply require that P(ξ 1 ∈ C) > 0, but our third condition also excludes the uninteresting cases where {S 1 , S 2 , . . . , S n ∈ C} = {S 1 = S 2 = • • • = S n = 0} almost surely. We consider the first exit time of the random walk from the cone defined by

T C = inf{n ≥ 1 : S n / ∈ C},
and wish to investigate the asymptotic distribution of (S 1 , . . . , S n ) conditional on {T C > n} as n → ∞.

We denote by C 1 the space of all continuous functions w : [0, 1] → R d , endowed with the topology of the uniform convergence and the corresponding Borel σ-algebra. Weak convergence of probability measures on C 1 will be denoted by the symbol ⇒.

Let Q n denote the distribution on C 1 of the process S n conditional on

{T C > n}, that is, for any Borel set B of C 1 , Q n (B) = P(S n ∈ B|T C > n). Note that, since C is a convex cone, this is equivalent to conditioning S n on {τ C (S n ) > 1}, where τ C (w) = inf{t > 0 : w(t) / ∈ C}, w ∈ C 1 .
We are interested in the weak convergence of the sequence of conditional distributions (Q n ). The one-dimensional case, where C = [0, ∞), has been investigated in the 60's and the 70's by many authors. It was Spitzer [START_REF] Spitzer | A Tauberian theorem and its probability interpretation[END_REF] who first announced a central limit theorem for the random walk conditioned to stay positive:

Q n (w(1) ≤ x) → 1 -exp(-x 2 /2), x ≥ 0.
But, apparently, he never published the proof. Note that the limit is the Rayleigh distribution. A first proof of the weak convergence of Q n was given by Iglehart in [START_REF] Iglehart | Functional central limit theorems for random walks conditioned to stay positive[END_REF] under the assumptions E(|ξ i | 3 ) < ∞ and ξ i nonlattice or integer valued with span 1. The limit is found to be the distribution of Brownian meander. Then Bolthausen proved in [START_REF] Bolthausen | On a functional central limit theorem for random walks conditioned to stay positive[END_REF] that these extra assumptions were superfluous. For the reader who may not be familiar with the Brownian meander, we will use a theorem of Durrett, Iglehart and Miller [START_REF] Durrett | Weak convergence to Brownian meander and Brownian excursion[END_REF] as a definition. Let W x be the distribution of the standard Brownian motion started at x. For any x > 0, we denote by

M x the distribution W x conditional on {τ C > 1}, that is M x (B) = W x (B|τ C > 1)
for any Borel set B of C 1 . Here, the definition of conditional probabilities is elementary since W x (τ C > 1) is positive for all x > 0. The distribution M of Brownian meander is the weak limit of M x as x → 0 + (see [START_REF] Durrett | Weak convergence to Brownian meander and Brownian excursion[END_REF], Theorem 2.1). Note that the existence of a limit is not straightforward since W 0 (τ C > 1) = 0. But in a sense, the Brownian meander is a Brownian motion started at 0 and conditioned to stay positive for a unit of time. The Brownian meander can alternatively be obtained by some path transformations of Brownian motion. Namely, it is the first positive excursion of Brownian motion with a lifetime greater than 1; it is also the absolute value of the rescaled section of Brownian motion observed on the interval [h, 1], where h is its last zero before t = 1 (see [START_REF] Bolthausen | On a functional central limit theorem for random walks conditioned to stay positive[END_REF] and [START_REF] Durrett | Weak convergence to Brownian meander and Brownian excursion[END_REF]).

With this in mind, the weak convergence of Q n to M can be stated in the following imprecise but intuitive way: the random walk conditioned to stay positive converges to a Brownian motion conditioned to stay positive.

We now turn to the two-dimensional case. If Q n does converge weakly, then its limit should naturally be the distribution of a Brownian motion conditioned to stay in the cone C for a unit of time. Such a process can be defined as the weak limit of conditioned Brownian motion in the same way as Brownian meander. As above, for x ∈ C o , let M x be the distribution of Brownian motion started at x and conditioned to stay in C for a unit of time. The following theorem is due to Shimura [START_REF] Shimura | Excursions in a cone for two-dimensional Brownian motion[END_REF] and has been extended in [START_REF] Garbit | Brownian motion conditioned to stay in a cone[END_REF] to any dimension d ≥ 2 for smooth cones.

Theorem 1.1 ([9], Theorem 2). As x ∈ C o → 0, the distribution M x converges weakly to a limit M .

The limit distribution M in this theorem will be referred to as the distribution of the Brownian meander (of the cone C). We will give more details about M in Section 2.

We now come to the main result of the present paper. We recall that a sequence (u n ) of positive numbers is regularly varying if it can be written as u n = n -α l n , where α ∈ R and (l n ) is slowly varying, i.e. lim n l [nt] /l n = 1 for all t > 0 (see for example [START_REF] Bojanic | A unified theory of regularly varying sequences[END_REF]). The exponent α is unique and called the index of regular variation. A non-increasing sequence of positive numbers (u n ) will be called dominatedly varying 1 if lim sup n u [nt] /u n is finite for all t ∈ (0, 1]. Theorem 1.2. Assume that the two-dimensional random walk has bounded increments. Then, the sequence of conditional distributions (Q n ) converges weakly to the Brownian meander if and only if P(T C > n) is dominatedly 1 This is strictly weaker than regular variation. For example, since n (1 + 1/n) is divergent, it is possible to construct a sequence of numbers 1 ≤ cn ≤ 2 such that :(i) for all n ≥ 1, cn+1 ≤ (1 + 1/n)cn, and (ii) for all ǫ > 0, there exist infinitely many n such that cn ≥ 2 -ǫ and cn+1 = 1. Then, the sequence un = cn/n is non-increasing and dominatedly varying, but not regularly varying since lim inf un+1/un ≤ 1/2 is not equal to 1 as it should be (see [START_REF] Weissman | A note on Bojanic-Seneta theory of regularly varying sequences[END_REF] for a very nice proof of this).

varying. In that case, P(T C > n) is regularly varying with index π/(2β), where β is the angle of the cone.

The assumption of bounded increments is only used in the proof of the tightness of (Q n ) which is taken from the paper [START_REF] Shimura | A limit theorem for two-dimensional random walk conditioned to stay in a cone[END_REF] of Shimura. We will discuss some extensions to the case where the increments are not bounded in Section 3.1. However, the rest of the proof of Theorem 1.2, which consists in a study of the (eventual) limit points of the sequence (Q n ), is completely independant of the assumption of bounded increments. Thus, we could have stated a more general (but not very useful) theorem by simply assuming that (Q n ) is tight. In order to avoid any confusion, the reader is advised that in any of the lemmas, propositions or theorems of this paper, the random walk (S n ) is not assumed to have bounded increments unless it is written explicitly.

Our Theorem 1.2 can be regarded as an extension of a previous result due to Shimura ([10], Theorem 1). Indeed, he proved that Q n ⇒ M if the distribution of the increments satisfies the following condition: there exists an orthogonal basis { u, v} of R 2 with v ∈ C o such that E(V |U ) = 0, where (U, V ) denotes the coordinates of ξ 1 in the new basis. But this condition does not seem to be very natural. For example, consider the simple random walk on Z 2 . It is well known that the coordinates (U, V ) of ξ 1 in the basis { u, v} = {(1, -1), (1, 1)} are independent, therefore Shimura's theorem applies to any of the cones {(x, y) : 0 ≤ y ≤ rx}, r > 1. However, if C is the octant {(x, y) : 0 ≤ y ≤ x}, there is no v in C o which satisfies the assumption of his theorem.

By comparison, regarding the example of the simple random walk, our Theorem 1.2 combined with the precise estimates of P(T C > n) given by Varopoulos in [START_REF] Varopoulos | Potential theory in conical domains[END_REF] shows that the weak convergence Q n ⇒ M holds for every adapted cone (Theorem 1.4). Indeed, Varopoulos estimates enable us to state an invariance principle (Theorem 1.3) that holds for a large class of random walks and adapted cones.

1.2.

Examples. Let L be the set of probability measures µ on R 2 with bounded support, mean zero, covariance matrix σ 2 I 2 with σ 2 > 0, and such that either µ has its support on Z 2 , or µ is absolutely continuous with respect to Lebesgue measure. Assume that the distribution µ of ξ 1 belongs to L. Recall S = (S n ) n≥1 is the associated random walk and let suppS denote its support.

If suppS ⊂ Z 2 , define p(n, x, y) by

p(n, x, y) = P x (S n = y; T C > n), x, y ∈ Z 2 ,
where P x (S ∈ * ) stands for P(x + S ∈ * ) as usual.

If µ is absolutely continuous, define p(n, x, y) by

p(n, x, y)dy = P x (S n ∈ dy; T C > n), x, y ∈ R 2 ,
that is y → p(n, x, y) is the density of the measure P x (S n ∈ * ;

T C > n).
Following the terminology of [START_REF] Varopoulos | Potential theory in conical domains[END_REF], we shall say that an adapted cone C is in general position with respect to µ if A) for all a > 0, there exist n ≥ 1 and ǫ > 0 such that ∀x, y ∈ C ∩ suppS, xy ≤ a implies p(n, x, y) ≥ ǫ.

B) for all a > 0, there exist n ≥ 1 and ǫ > 0 such that

∀x ∈ C ∩ suppS, d(x, ∂C) ≤ a implies P x (T C ≤ n) ≥ ǫ.
Assuming that C is in general position with respect to µ, Varopoulos obtained in [START_REF] Varopoulos | Potential theory in conical domains[END_REF] precise estimates of the tail distribution of T C . In particular, his results show that there exists γ > 1 such that

γ -1 n -π/2β ≤ P(T C > n) ≤ γn -π/2β ,
where β is the angle of the cone. This, combined with Theorem 1.2, proves the following invariance principle.

Theorem 1.3. Assume that the distribution of ξ 1 belongs to L, and that C is in general position. Then (Q n ) converges weakly to M .

The reader may have noticed that the simple random walk on Z 2 does never fit the condition A) above because of a parity problem: if y and z are two neighbours, then p(n, x, y) and p(n, x, z) can not be positive simultaneously. However, it is possible to get around this and find a good estimate of P(T C > n) in this context. Let us explain how.

Let (S n ) be the simple random walk on Z 2 and let C be an adapted cone. For any x ∈ Z 2 ∩ C, one can always find a path (on Z 2 ) that joins 0 to x without leaving C, hence there exists k ≥ 0 such that P(S k = x, T C > k) > 0. Fix x ∈ Z 2 deep inside the cone so that d(x + C, ∂C) > 1. Since the steps of (S n ) are bounded by 1, we have {S 2i ∈ C for all i = 1 . . . m} ⊂ {S j ∈ Cx for all j = 1 . . . 2m}, and consequently

P(S 2i ∈ C for all i = 1 . . . m) ≤ P x (S j ∈ C for all j = 1 . . . 2m). Now, choose k such that α := P(S k = x, T C > k) > 0 and set m = [ n-k 2 ] + 1.
Using the Markov property of the random walk and the preceding inequality, we obtain, for n > k,

P(T C > n) ≥ P(S k = x, T C > n) ≥ P(S k = x, T C > k)P x (S 1 , S 2 , . . . , S 2m ∈ C) ≥ αP(S 2 , S 4 , . . . , S 2m ∈ C). (1)
On the other hand, we clearly have ( 2)

P(T C > n) ≤ P(S 2 , S 4 , . . . , S 2m-2 ∈ C).
Hence, the problem reduces to that of finding a good estimate of P(S 2 , S 4 , . . . , S 2m ∈ C); the benefit of this reduction is that the cone C is in general position with respect to the random walk (S 2m ) m≥0 . One of the key point is that the increments of (S 2m ) m≥0 are equal to 0 with positive probability so that the parity problem does not occur. The complete proof is not difficult, so we omit it. Now, the combination of Varopoulos estimate for P(S 2 , S 4 , . . . , S 2m ∈ C) and inequalities ( 1) and ( 2) gives the expected result :

γ -1 n -π/2β ≤ P(T C > n) ≤ γn -π/2β
for some γ > 1. Therefore, by application of our Theorem 1.2 we obtain :

Theorem 1.4. Let (S n ) be the simple random walk on Z 2 and C be an adapted cone. Then Q n ⇒ M .

1.3. Degenerated cases. The question rises whether the sequence of conditional distributions (Q n ) can converge to some limit Q = M . In that case P(T C > n) would not be dominatedly varying. We do not know any example in dimension 2 (with C adapted to the random walk) but there are some in dimension 3.

Example 1.5. Let (S n ) be the simple random walk on Z 3 , and take

C = {(x, y, z) ∈ R 3 : 0 ≤ x/2 ≤ y ≤ 2x}. Conditional on {T C > n}, the process {S k , k = 1 . . . n} is a simple random walk on the axis {(0, 0, z) : z ∈ Z}. Thus Q n ⇒ Q, where Q is the law of the process {(0, 0, √ 3B t ), t ≥ 0}, B t being a stantard one-dimensional Brownian motion. Here P(T C > n) = (1/3) n .
Example 1.6. Here again consider the simple random walk on Z 3 , and take C = {(x, y, z) ∈ R 3 : 0 ≤ x/2 ≤ y ≤ 2x and z ≥ 0}. Conditional on {T C > n}, the process {S k , k = 1 . . . n} is a simple random walk on the axis {(0, 0, z) : z ∈ Z} conditioned to stay positive. Hence Q n ⇒ Q, where Q is the law of the process {(0, 0, √ 3M t ), t ≥ 0}, M t being a one-dimensional Brownian meander. Here

P(T C > n) ∼ (1/3) n (πn) -1/2 .
These examples show that it is possible that Q n converges to some limit even if P(T C > n) is not dominatedly varying. But then, the limit process is degenerated. The following proposition, which holds in any dimension and for any adapted cone, shows that it is a general fact.

Proposition 1.7. Suppose Q n ⇒ Q weakly on C 1 . If P(T C > n)
is not dominatedly varying then the limit process lives on the boundary of C, i.e.

Q(∀t ∈ [0, 1], w(t) ∈ ∂C) = 1.
The proofs of Theorem 1.2 and Proposition 1.7 are deferred to Section 3. Prerequisites are collected in Section 2.

Preparatory material

We collect here some of the results we need to prove Theorem 1.2.

2.1. More on Brownian meander. Let (X t , 0 ≤ t ≤ 1) be the canonical process on C 1 for which X t (w) = w(t).

Given 0 < β ≤ π, let C ⊂ R 2 be the cone {(r cos θ, r sin θ) : r > 0 and 0 < θ < β}. (One can replace the condition 0 < θ < β in the definition of C by 0 ≤ θ < β, 0 < θ ≤ β or 0 ≤ θ ≤ β, for the exit times of these cones are almost surely equal relative to Wiener distribution; hence the distribution of Brownian motion conditional on {τ C > 1} does not depend on that choice.) Let us denote (as before) M the distribution on C 1 of the Brownian meander of the cone C. The Brownian meander is a continuous, non-homogeneous Markov process, with transition density given by M (X t ∈ dy) = e(t, y)dy

(3) = r 2α 2 α Γ(α)t 2α+1 exp(- r 2 2t
) sin(2αθ)W y (τ C > 1t)dy for 0 < t ≤ 1 and y = (r cos θ, r sin θ), r > 0, 0 ≤ θ ≤ β; α = π/(2β). And M (X t ∈ dy|X s = x) = p(s, x, t, y)dy

= p C (t -s, x, y) W y (τ C > 1 -t) W x (τ C > 1 -s) dy
for 0 < s < t ≤ 1 and x, y ∈ C o , where p C (t, x, y) is the heat kernel of C.

(See [START_REF] Garbit | Brownian motion conditioned to stay in a cone[END_REF] for a derivation of these formulas; the formula (3.2) given in [START_REF] Shimura | Excursions in a cone for two-dimensional Brownian motion[END_REF] for the transition density was misprinted.) From formula (3), one can already obtain a necessary condition for the weak convergence of the conditional distribution Q n to the Brownian meander. The conclusion of the following proposition was given by Shimura in [10] as a consequence of his limit theorem.

Proposition 2.1. If Q n ⇒ M in C 1 , then P(T C > n) is regularly varying with index α = π/(2β). Proof. Given t > 1, define φ n (x) = P x √ n (T C > [nt] -n).
By the Markov property of the random walk and the definition of Q n , we have

P(T C > [nt]|T C > n) = Q n (φ n (X 1 )).
If x ∈ C o and x n → x, then it follows from Donsker's theorem and the Portmanteau theorem ([1], Theorem 2.1) that

φ n (x n ) → φ(x) = W x (τ C > t -1).
Since Q n ⇒ M , and X 1 ∈ C o M -a.s., the continuous mapping theorem ([1], Theorem 5.5) shows that

Q n (φ n (X 1 )) → M (φ(X 1 )).
Using (3), the limit can be expressed as

C e(1, x)W x (τ C > t -1)dx = C te(1, √ ty)W y (τ C > 1 -1/t)dy,
where we have made the change of variables x = √ ty and used the scaling invariance of Brownian motion. But the last integrand is equal to t -α e(1/t, y), therefore M (φ(X 1 )) = t -α . This proves that ( 4)

P(T C > [nt])/P(T C > n) → t -α
for all t > 1. In a similar way, it can be proved that (4) also holds for all t ∈ (0, 1), thus P(T C > n) is regularly varying with index α = π/(2β).

2.2. More on conditioned random walk. Let us finally state here without proof two easy but important facts about the conditioned random walk. The first one is the Markov property which is inherited from the original unconditioned random walk. Let F t be the σ-algebra generated by the random variables {X s , s ≤ t}. The shift operator θ t on C 1 is defined by θ t (w)(s) = w(t + s).

For all x ∈ C and t ∈ [0, 1], we denote by Q x,t n the distribution of x + S n conditional on {τ C (x + S n ) > t}. Then we have:

Proposition 2.2. Let t = k/n ∈ [0, 1] be given. For any A ∈ F t and B ∈ C 1 , Q n (A; θ -1 t B) = Q n (A; Q Xt,1-t n (B)).
The second fact is a limit theorem for the conditioned normalized random walk started inside the cone. For x ∈ C o and t ∈ (0, 1], let M x,t denote the distribution of the standard Brownian motion started at x and conditioned to stay in C until time t, that is

M x,t (B) = W x (B|τ C > t),
for any Borel set B of C 1 . If t n = k n /n → t and x n → x, then it follows from Donsker's theorem and the Portmanteau theorem that

P(x n + S n ∈ B; τ C (x n + S n ) > t n ) → W x (B; τ C > t), for any Borel set B such that ∂B is W x -negligible. Since W x (τ C > t) > 0 (because x ∈ C o ), we obtain: Theorem 2.3. Let t ∈ (0, 1] and x ∈ C o . If t n = k n /n → t and x n → x, then Q xn,tn n ⇒ M x,t .
Detailed proofs of Proposition 2.2 and Theorem 2.3 can be found in [START_REF] Garbit | Contributions à l'étude d'une marche aléatoire centrifuge et théorèmes limites pour des processus aléatoires conditionnés[END_REF].

Proofs of results

The proof of Theorem 1.2 begins with the analysis of the approach developped by Shimura in [START_REF] Shimura | A limit theorem for two-dimensional random walk conditioned to stay in a cone[END_REF]: he proved that the sequence (Q n ) is tight and then observed (without stating it as a general fact) that a sufficient (and necessary) condition for the convergence Q n ⇒ M is that any weak limit point Q of (Q n ) satisfies Q(w(t) ∈ ∂C) = 0 for all t ∈ (0, 1]. We explain this in Section 3.1. We will then show in Section 3.2 how this condition relates to the tails of the exit time. This is the main novelty of this paper.

3.1. Shimura's approach. The paper [START_REF] Shimura | A limit theorem for two-dimensional random walk conditioned to stay in a cone[END_REF] of Shimura contains two significant results. The first one is the tightness of the sequence (Q n ) under the assumption that the increments of (S n ) are bounded. Shimura's proof of tightness is quite technical and requires an extension of Theorem 2.3 to sequences x n → x ∈ ∂C \ {0}. We do not recall it here and refer the reader to [START_REF] Shimura | A limit theorem for two-dimensional random walk conditioned to stay in a cone[END_REF], Lemma 3. Let us simply mention that the assumption of bounded increments can be replaced by the following one :

(5) P( max

i=1...n ξ i > √ n|T C > n) → 0.
This is explained in [START_REF] Garbit | Contributions à l'étude d'une marche aléatoire centrifuge et théorèmes limites pour des processus aléatoires conditionnés[END_REF]. Of course, condition ( 5) is satisfied if the steps ξ i are bounded, but in other cases it is not clear if one should impose an additional moment condition or not. For example, suppose that we already know that P(T C > n) ≥ γn -α for some positive constants γ and α 2 . A trivial upper bound is given by :

P( max i=1...n ξ i > √ n|T C > n) ≤ P(max i=1...n ξ i > √ n) P(T C > n) ≤ γ -1 n α+1 P( ξ 1 > √ n).
Hence condition ( 5) is satisfied whenever E( ξ 1 2α+2 ) is finite. Since α should not be less than π/2β where β is the angle of C, this condition on moments asks for a finite third moment in the half-plane case (for which we already know that a second moment is sufficient since it is similar to the one-dimensional case), a finite fourth moment in the quarter-plane case, or a finite sixth moment in the octant case. Clearly, we are far from an optimal condition! In the same spirit, if one knows a priori nothing about the asymptotic behavior of P(T C > n), it is still possible to obtain (5) under a very strong integrability assumption. Namely, the condition is

(6) E ξ 1 2 a ξ 1 2 < ∞,
where a = P(ξ 1 ∈ C) -1 . Indeed, by a classical argument, condition [START_REF] Garbit | Brownian motion conditioned to stay in a cone[END_REF] implies that na n P(

ξ 1 > √ n) → 0. But, since C is a semi-group, we also have P(T C > n) ≥ a -n . Therefore, P( max i=1...n ξ i > √ n|T C > n) ≤ nP( ξ 1 > √ n) P(T C > n) ≤ na n P( ξ 1 > √ n) → 0,
and consequently, the sequence (Q n ) is tight. We insist on the fact that all these considerations only concern the question of the tightness of (Q n ).

In the rest of this section, we do not assume that the random walk (S n ) has bounded increments, nor any k-th moment for k > 2.

2 Apart from the random walks with bounded increments that were considered in our Section 1.2, this lower bound is also known for a certain class of random walks with unbounded increments, see Section 6 of [START_REF] Varopoulos | Potential theory in conical domains[END_REF].

We now turn to the study of the eventual limit points of (Q n ). Let Q be a probability measure on C 1 . We shall say that ∂C is unreachable for Q if Q(w(t) ∈ ∂C) = 0 for all t ∈ (0, 1]. By virtue of (3), the boundary of C is unreachable for M . Perhaps, the most striking result in Shimura's paper is that M is the only possible subsequential limit of (Q n ) for which ∂C is unreachable. The proof is illuminating and we believe it is of interest to reproduce it here. Proposition 3.1 ([10], Proof of Theorem 1). Let Q be a limit point of the sequence (Q n ). Then Q = M if and only if ∂C is unreachable for Q.

Proof. Let Q be a weak limit point of the sequence (Q n ). There exists a subsequence (Q n ′ ) which converges weakly to Q. To simplify the notation, we shall suppose that the whole sequence (Q n ) converges to Q. We will show that Q and M have the same one-dimensional distributions; the generalization to other finite-dimensional distributions is straightforward.

Fix t ∈ (0, 1] and let f be a bounded continuous real function. We have to show that

Q(f (X t )) = M (f (X t )).
First choose 0 < λ < t and set λ n = [nλ]/n and t n = [nt]/n. Note that λ n → λ and t n → t as n → ∞. Given a vector u ∈ C o we set

C ǫ = ǫu + C and ∆ ǫ = C \ C o ǫ . Note that ∩ ǫ>0 ∆ ǫ = ∂C. For all ǫ > 0, define J n ǫ = Q n (X λn ∈ C ǫ ; f (X tn )). Then, |Q n (f (X tn )) -J n ǫ | ≤ KQ n (X λn ∈ ∆ ǫ )
, where K is a bound for |f |. By the continuous mapping theorem Q n (X λn ∈ dx) ⇒ Q(X λ ∈ dx). A standard use of the Portmanteau theorem then shows that lim

ǫ→0 lim sup n→∞ Q n (X λn ∈ ∆ ǫ ) ≤ Q(X λ ∈ ∂C) = 0 .
In addition, (continuous mapping theorem again) we have

lim n→∞ Q n (f (X tn )) = Q(f (X t )).
Therefore, lim

ǫ→0 lim sup n→∞ |Q(f (X t )) -J n ǫ | = 0.
Thus, it remains to prove that lim

ǫ→0 lim n→∞ J n ǫ = M (f (X t )).
By the Markov property of Q n , we have

J n ǫ = Q n X λn ∈ C ǫ ; Q X λn ,1-λn n (f (X tn-λn )) = Q n (X λn ∈ C ǫ ; φ n (X λn )), where φ n (x) = Q x,1-λn n (f (X tn-λn )), x ∈ C. If x n → x ∈ C o , then by Theorem 2.3, lim n→∞ φ n (x n ) = M x,1-λ (f (X t-λ )) =: φ λ (x).
Hence, if w is such that w(λ) / ∈ ∂C ǫ , and if w n → w uniformly on [0, 1], then [START_REF] Iglehart | Functional central limit theorems for random walks conditioned to stay positive[END_REF] lim

n→∞ 1 1 Cǫ (w n (λ n ))φ n (w n (λ n )) = 1 1 Cǫ (w(λ))φ(w(λ)).
Let S be the set of all ǫ > 0 such that Q(X λ ∈ ∂C ǫ ) = 0. The set (0, ∞) \ S is at most countable. Fix ǫ ∈ S. By the continuous mapping theorem, it follows from [START_REF] Iglehart | Functional central limit theorems for random walks conditioned to stay positive[END_REF] that

lim n→∞ J n ǫ = Q(X λ ∈ C ǫ ; φ λ (X λ )). Now, letting ǫ → 0 through S gives lim ǫ→0 lim n→∞ J n ǫ = Q(X λ ∈ C o ; φ λ (X λ )) = Q(φ λ (X λ )),
since Q(X λ ∈ ∂C) = 0. This last expression does not depend on λ, so we are going to let λ → 0. To do this, select a sequence λ n → 0. Then, with Qprobability one, w(λ n ) ∈ C o for all n (by hypothesis), and w(λ n ) → w(0) = 0 (by continuity of the paths). Therefore, an easy modification of Theorem 1.1 using only the scaling invariance of Brownian motion shows that φ λn (w(λ n )) = M w(λn),1-λn (f (X t-λn )) → M (f (X t )).

This holds for Q-almost all w, so by the dominated convergence theorem, we have lim

n→∞ Q(φ λn (X λn )) = M (f (X t )),
which completes the proof.

Remark 3.2. We point out that Proposition 3.1 holds in any dimension provided that the Brownian meander, defined as the weak limit of Brownian motion conditioned to stay in the cone for a unit of time, exists. For example, this holds for cones with a smooth boundary (see [START_REF] Garbit | Brownian motion conditioned to stay in a cone[END_REF] for further details).

Tails of the exit time and unreachability of the boundary.

In this section, the dimension is an arbitrary integer d ≥ 1, and C is an adapted cone of R d . We do not assume that the random walk (S n ) has bounded increments. It turns out that the unreachability of the boundary of C is closely related to the asymptotic behavior of the tail distribution of the exit time T C . The first lemma gives a sufficient condition for unreachability.

Lemma 3.3. Let Q be a limit point of the sequence (Q n ). If P(T C > n) is dominatedly varying, then ∂C is unreachable for Q.

Proof. Without loss of generality, we shall assume that σ 2 = 1. Choose a vector u ∈ C o and set C ǫ = ǫu+C and ∆

ǫ = C \C o ǫ . Note that ∩ ǫ>0 ∆ ǫ = ∂C. Fix 0 < s < t ≤ 1 and define p(n, ǫ, R) = P( S [ns] ≤ √ nR; S [nt] ∈ √ n∆ ǫ |T C > n), n, ǫ, R > 0.
We shall first prove that By the Markov property of the random walk, we have

P(T C > n; S [ns] ≤ √ nR; S [nt] ∈ √ n∆ ǫ ) ≤ P(T C > [ns]; S [ns] ≤ √ nR; P S [ns] (S kn ∈ √ n∆ ǫ )) ≤ P(T C > [ns]) sup{P x (S kn ∈ √ n∆ ǫ )) : x ≤ √ nR}
where

k n = [nt] -[ns]
. Thus, with our choice of γ, we obtain

(9) p(n, ǫ, R) ≤ γ sup{P(z + S kn / √ n ∈ ∆ ǫ )) : z ≤ R}. If z n → z,
Q n ( w(s n ) ≤ R; w(t n ) ∈ ∆ ǫ ) = 0.
By the Portmanteau theorem, this implies

∀R > 0, Q( w(s) ≤ R; w(t) ∈ ∂C) = 0.
Letting R → ∞ completes the proof.

If the whole sequence (Q n ) converges weakly to some limit Q, a converse to Lemma 3.3 holds: Lemma 3.4. Suppose Q n ⇒ Q in C 1 . If P(T C > n) is not dominatedly varying, then Q(w(t) ∈ ∂C) = 1 for all t ∈ (0, 1].

Proof. Here again, we shall assume without loss of generality that σ 2 = 1. Suppose there exists t ∈ (0, 1] such that Q(w(t) ∈ ∂C) < 1 and set t n = [nt]/n. Define C ǫ and ∆ ǫ as in the proof of Lemma 3.3. Then, by the Portmanteau theorem, lim sup

n→∞ Q n (w(t n ) ∈ ∆ ǫ ) ≤ Q(w(t) ∈ ∆ ǫ ).
Hence, there exists ǫ > 0 such that [START_REF] Shimura | A limit theorem for two-dimensional random walk conditioned to stay in a cone[END_REF] lim sup

n→∞ Q n (w(t n ) ∈ ∆ ǫ ) < 1.
We shall prove that this implies the dominated variation of P(T C > n). Proof of Theorem 1.2. Suppose that the random walk (S n ) has bounded increments. Then, by Lemma 3 in [START_REF] Shimura | A limit theorem for two-dimensional random walk conditioned to stay in a cone[END_REF], we know that the sequence (Q n ) is tight on C 1 . If P(T C > n) is dominatedly varying, then it follows from Theorem 3.1 and Lemma 3.3 that M is the only possible limit point of (Q n ). Therefore Q n ⇒ M . The converse part is Proposition 2.1, which also gives the index of regular variation α = π/2β. Note that the dominated variation of P(T C > n) can be derived from Lemma 3.4, since we know from (3) that M (w(t) ∈ ∂C) = 0 for all t ∈ (0, 1].

Proof of Proposition 1.7. Assume the hypotheses of Proposition 1.7 are satisfied. Then Lemma 3.4 ensures that Q(w(t) ∈ ∂C) = 1 for all t ∈ (0, 1]. Thus Q(∀t ∈ Q ∩ [0, 1], w(t) ∈ ∂C) = 1. The result follows since every path w ∈ C 1 is continuous, Q is dense in [0, 1] and ∂C is a closed set.

  , ǫ, R) = 0.Since P(T C > n) is dominatedly varying, there exists a positive constant γ such that ∀n ≥ 1,P(T C > [ns]) ≤ γP(T C > n).

3 . 3 .

 33 Let s ∈ (0, 1) be given and set m =[ns]. By the Markov property of the random walk, we haveP(T C > n) ≥ P(T C > n; S [mt] ∈ √ mC ǫ ) ≥ P(T C > [mt]; S [mt] ∈ √ mC ǫ ; P S [mt] (T C > n -[mt])) ≥ P(T C > m; S [mt] ∈ √ mC ǫ ) inf x∈ √ mCǫ P x (T C > n -[mt]). (11) But, since C is a semi-group, p m := inf x∈ √ mCǫ P x (T C > n -[mt]) ≥ P √ mǫu (T C > n -[mt]),and it follows from Donsker's invariance principle and the Portmanteau theorem that(12) lim inf n→∞ p m ≥ p = W ǫu (τ C o > s -1t) > 0.Now, dividing both sides of (11) by P(T C > m) givesP(T C > n) P(T C > m) ≥ Q m (w(t m ) ∈ C ǫ ) × p m .Thus, by virtue of (10) and[START_REF] Varopoulos | Potential theory in conical domains[END_REF], we obtainlim inf n→∞ P(T C > n) P(T C > m) > 0.Therefore P(T C > n) is dominatedly varying. Final steps. Let us finally give the proofs of our main results, namely Theorem 1.2 and Proposition 1.7.

Acknowledgments

I would like to thank E. Lesigne, M. Peigné and an anonymous referee for helping in many ways to improve the original manuscript.

This work was partially supported by Agence Nationale de la Recherche Grant ANR-09-BLAN-0084601.


J Aubry

A H Ben

Ahmed Is With Satie

B Multon

Satie

Keywords: Bi-objective optimization, design methodology, direct-drive, operating profiles, optimization methods, particle swarm, permanent-magnet generators, sizing optimization, wave energy conversion. I

This article sets forth a design optimization methodology for an SMPM machine drive on an operating profile. The optimization parameters serve to describe both the machine geometry and electrical ratings of the electronic power converter. Each operating point is treated independently, and current control is optimized at every operating point in order to not only minimize machine drive losses but also satisfy a number of constraints. This optimization methodology will then be applied to the design of a direct-drive conversion chain for a wave energy converter (WEC).

INTRODUCTION

EVERAL approaches to designing a permanent magnet machine that make use of optimization algorithms, either deterministic [START_REF] Slemon | Modeling and Design Optimization of Permanent Magnet Motors[END_REF][START_REF] Messine | Optimal Design of Electromechanical Actuators: A New Method Based on Global Optimization[END_REF] or stochastic [START_REF] Bianchi | Design optimization of electric motor by genetic algorithms[END_REF][START_REF] Boules | Design Optimization of Permanent Magnet DC Motors[END_REF], have been presented in the literature. Very little research has focused however on such a design exercise over a given operating profile [START_REF] Biais | Optimization of a permanent magnet traction motor for electric vehicle[END_REF]. Yet, for an increasing number of applications (electric traction, renewable energy production, etc.), it proves impossible to define a typical operating point for a given design.

The present article proposes developing a design optimization methodology on the operating profile of an SMPM machine drive for a wave energy conversion application. The term "operating profile" refers herein to a series of torque-speed operating points, denoted , Ω , that the converter-machine set must achieve. All points on the operating profile are taken into account in order to determine objectives and satisfy constraints. The two competing optimization objectives selected in this case are the cost of the electric chain (machine active parts + power electronics) and electrical energy (generated in our case). Both the incorporation and optimization of electronic power converter ratings require adopting a current control strategy (direct and transverse axis current) for each operating point

, Ω .

This study will examine two contributions of flux weakening for permanent magnet machine drives:

-Extension of the constant power speed range. In several articles [START_REF] Jahns | Flux-Weakening Regime Operation of an Interior Permanent-Magnet Synchronous Motor Drive[END_REF][START_REF] Sneyers | Field Weakening in buried Permanent Magnet AC Motor Drives[END_REF][START_REF] Soong | Field-weakening performance of brushless synchronous AC motor drives[END_REF][START_REF] Morimoto | Expansion of operating limits for permanent magnet motor by current vector control considering inverter capacity[END_REF][START_REF] Sebastian | Operating Limits of Inverter-Driven Permanent Magnet Motor Drives[END_REF][START_REF] Multon | Analyse des possibilités de fonctionnement en régime de désexcitation des moteurs à aimants permanents[END_REF], it is shown that electronic flux weakening (through modification of the leading angle or injection of a demagnetizing direct axis current) makes it possible, for given converter ratings, to extend the operating speed range of a machine drive with respect to a conventional control (e.g. an 0 control).

-Minimization of losses [START_REF] Mademlis | Loss minimization in surface permanent magnet synchronous motors drives[END_REF][START_REF] Fernandez-Bernal | Model-based loss minimization for DC and AC vector-controlled motors including core saturation[END_REF][START_REF] Morimoto | Loss minimization control of permanent magnet synchronous motors drives[END_REF]. By restricting losses to just copper losses, then loss minimization at each operating point necessitates a so-called "maximum torque per ampere" control, which within the scope of non-salient machines, may be expressed very simply as 0. By adding other losses (magnetic, converter) to this formulation, the control law for minimizing losses at each point becomes more complex, while still remaining feasible.

A broad array of research has been conducted to take advantage of these two contributions from flux weakening. Yet such research has been limited to optimizing only the control strategy [START_REF] Chédot | Comparison of direct and adaptative optimal controls for interior permanent magnet synchronous integrated starter generator[END_REF], or optimizing the control and machine geometry (excluding the converter) [START_REF] De La Barrière | PM Axial Flux Machine Design for hybrid Traction[END_REF], or lastly optimizing the control-converter-machine set while focusing on just a single criterion (either the machine mass [START_REF] Kreuawan | Modeling and optimal design for railway applications[END_REF] or losses on the operating profile [START_REF] Figueroa Barnier | Modélisation des entraînements à grande plage de vitesse en vue de leur conception[END_REF]).

The first part of this paper will discuss the general optimization problem associated with the design of a machine drive; this will entail presenting the two selected competing objectives in addition to the set of machine drive design parameters and relevant constraints. The second part will demonstrate how to optimize the direct axis current in order to minimize a loss criterion for each operating point, with respect to both converter design parameters and magnet demagnetization. Afterwards, the models used to perform this optimization routine will be described in detail. The last part will apply the approach developed to the study of a pendulum system devoted to wave energy conversion.

II. PRESENTATION OF THE OPTIMIZATION PROBLEM

A. Optimization objectives From a classical perspective, the design objectives for a machine drive could consist of machine mass, machine cost or total losses. Two specific and competing objectives will be discussed here: total cost of the machine drive, and electrical energy.

1) Objective 1: Cost of the electric chain In order to evaluate the machine drive cost, we consider the sum of the power electronics cost (inverter, filters, control circuits) and machine cost, e.g.:

(1)

Our assumption is that the cost of the converter is a function of its apparent power (2) (stemming from observation). As for the machine cost, only the sum of costs for raw materials and active parts have been taken into account [START_REF] Bianchi | Design optimization of electric motor by genetic algorithms[END_REF], with numerical values being listed in Table II

Coefficients and may be adjusted depending on the application (large or small series), and they remain valid for a given power range and conversion structure. The coefficients for an IGBT three-phase bridge will be provided in the application example further below.

2) Objective 2: Electrical energy Depending on the specific application, e.g. whether motor or generator, the electrical energy may be either consumed or generated. If we adopt a motor convention, in both cases this electrical energy requires minimization:

Ω , (4)
where , are the total machine drive losses at operating point , Ω with a duration of .

B. Optimization parameters

The target to be optimized is an SMPM machine drive. The optimization parameters thus need to describe the two physical components making up this system : machine and converter.

Let's assume the DC bus voltage to be constant. The converter will then be parameterized by its rated current ( with 3 and /2√2). In order to effectively describe the machine, let's identify just the geometric parameters shown in Fig. 1 andlisted

C. Optimization constraints

To complete the description of this optimization problem, several constraints still need to be defined. Some constraints must be satisfied at each operating point , Ω :

-Rated voltage:

-Rated current:

-Saturation: 1,5

-Demagnetization: Other geometrical constraints however must be satisfied for each machine drive, e.g.:

-

0.2 3 - -
The thermal constraint must be treated with special attention. The temperature rise relative to ambient temperature must be calculated strictly by means of a transient thermal model. The maximum value over the entire increasing temperature profile must therefore remain below a threshold value :

ΔΘ ΔΘ Should the sequencing of operating points not be known, only one thermal model in the steady-state mode may be used and the temperature rise for each operating point is denoted ΔΘ . Two scenarios then become possible:

-When thermal time constants of the machine are large relative to the loss variation periods, the temperature rise constraint can be calculated based on average loss value : ΔΘ ΔΘ -When the thermal time constants are small relative to the loss variation periods, temperature rise may be considered as instantaneous and the corresponding constraint appears on the least favorable operating point :

ΔΘ ΔΘ Depending on the application and level of power under consideration, one of these two scenarios may be envisaged.

Otherwise, it would be necessary to calculate the temperature rise profile on the basis of a transient thermal model, which could cause a higher computation time.

III. CURRENT CONTROL OPTIMIZATION ON AN OPERATING POINT

In order to evaluate machine drive losses , at operating point , Ω , the direct axis current , still needs to be determined. The transverse axis current , , in case of a non-salient machine, is defined by the torque to be supplied.

The values of currents , at each operating point should constitute strictly optimization parameters to the same extent as machine drive design parameters. On an operating profile containing a large number of points however, this rigorous approach is infeasible since it leads to an excessive increase in the number of optimization parameters.

The alternative approach adopted herein consists of locally optimizing the value of current , for each operating point. Let's assign the minimization of total machine drive losses as the optimization criterion, which in turn contributes to optimizing one of the two design objectives, namely electrical energy . As indicated in the previous section however, some constraints must be satisfied at each operating point

, Ω : voltage, current, saturation and demagnetization. The value , must therefore be optimized subject to these constraints. Should no optimization solution exist, e.g. no , value found to satisfy all four constraints, then the particular operating point cannot be reached. An illustration of this optimization protocol is presented in Fig. 2

(only voltage and current constraints are illustrated).

A new constraint will now be defined that encompasses the four mentioned above: each operating point , Ω must be reached.

In practice, this constrained optimization routine is performed by systematically computing several discrete values lying within the interval ; 0 . The selected , value is the one that minimizes losses while satisfying constraints. The number of computed discrete values will be established on the basis of a compromise between level of accuracy and computation time.

IV. MACHINE DRIVE MODELIZATION

The goal of this modeling effort is to determine the objectives to be optimized based on a set of optimization parameters. The "economical cost" objective does not require any kind of physical model, but the electrical energy calculation requires knowing the losses at each operating point.

This modelization is restricted to the first harmonic and decomposes airgap induction on two axes that rotate synchronously with the electrical angular velocity ω (d-q axes). The hypotheses adopted to build this model consist of: infinite iron permeability, and the independence of electrical quantities with respect to temperature.

A. Electrical machine model

The equivalent circuits in both axes are depicted in Fig. 3. The iron loss resistance is neglected during the calculation of electrical values.

The decomposition into d-q axis has been chosen so that the RMS current value per phase and the RMS phase to neutral voltage are both expressed as follows:

1 √3

(5)

1 √3 (6)
Fig. 3: d-and q-equivalent circuits of the SMPM machine

Rs:

Armature winding resistance per phase Rf:

Equivalent iron loss resistance Φf: √3Φ Φe:

Flux linkage of the PM per phase (r.m.s. value) Lm:

Magnetizing inductance Ll:

Leakage inductance

The electromagnetic torque equals: Φ

B. Model of machine drive losses

The resistance of one phase is determined from the set of geometric parameters. Copper resistivity is set at the maximum allowable temperature. The dependence of resistivity or other physical quantities on temperature could be incorporated by means of an iterative calculation, yet computation time would rise excessively.

, , , (8)
Iron losses are calculated in the first harmonic hypothesis. The consideration of harmonic losses in the model would not cause any additional difficulties. Equation ( 9) is used to calculate loss densities:

, , , 2 |Ω | Ω 2 (9)
Total iron losses are derived from ( 9) by summing total losses in both the teeth and stator yoke. Loss densities are considered to be homogeneous in each of these two parts (rotor losses are assumed to equal zero for this model in the first harmonic under steady-state mode).

Losses in the electronic power converter are composed of the sum of conduction losses [START_REF] Sebastian | Operating Limits of Inverter-Driven Permanent Magnet Motor Drives[END_REF] and switching losses (13) within both the diodes and transistors (in this case IGBT).

, , , ,

where x is the component under consideration (diode or IGBT), , the threshold voltage of the component and , its dynamic resistance. Observations gleaned from manufacturers' documentation enable deducing, for a given maximum voltage rating, the scale law as a function of the maximum current rating of the various components:

, / √2 (11)
, / √2

Switching losses are assumed to be proportional to the switched current:

√2/ (13)
where and are the switching loss coefficients in the IGBT and diode, respectively. The hypothesis will be adopted whereby for a given maximum voltage rating, these coefficients remain independent of the maximum current rating of the switches.

C. Magnetic machine model

An expression of the magnetizing inductance is shown in ( 14) for a diametral winding with one slot per pole and per phase:

3 2 4 ( 14
)
where is Carter's coefficient [START_REF] Gieras | Permanent magnet motor technology: Design and application[END_REF]. The leakage inductance is obtained from a calculation of the magnetic energy stored in a single slot (end windings and tooth tips have been neglected herein), which yields:

2 3 (15)
D. Thermal machine model The equivalent thermal circuit of the machine is displayed in Fig. 4. Due to symmetry considerations, only half of a tooth pitch has been examined.

To begin, the curvature effects are neglected and losses are distributed; in this case, the calculation of equivalent circuit elements requires special attention [START_REF] Gerling | Novel lumped-parameter thermal model for electrical systems[END_REF]. The expression for thermal resistances is detailed for the slot only as follows:

_ , 1 2 (16) _ , 1 2 (17)
, corresponds to the copper loss produced inside a half-slot.

corresponds to the equivalent thermal conductivity of the stator winding (copper + insulation), set equal to 0,5 .

. in the case of vacuum impregnation. The resistances , and , enable incorporating the tangential heat flux from the slot to the tooth; this flux cannot be neglected in the case where slot depth is greater than slot width.

The resistance serves to model the thermal exchange taking place by convection and radiation.

The value of the temperature rise constraint ΔΘ is calculated inside the winding; it corresponds to the temperature gain between the hottest point of the slot and the ambient temperature.

For the sake of simplification, a thermal model of the electronic power converter was not developed since maximum current ratings have been established on the basis of manufacturer guidelines; moreover, it can be considered that cooling will be adapted in order to maintain the junction temperature at an acceptable value.

V. DESIGN OPTIMIZATION ON AN OPERATING PROFILE

This article has demonstrated how to evaluate and optimize the performance of a given machine drive at an operating point. Fig. 5 offers the optimization flowchart. The cost of the electric chain, along with the geometric constraints, may be directly assessed from the parameters presented. Electrical energy and the thermal constraint are calculated from losses. If for a given operating point no solution exists that optimizes , , then the corresponding set of optimization parameters will not be selected as part of the overall design optimization.

A. Bi-objective optimization algorithm

We introduced a bi-objective optimization algorithm, based on both Particle Swarm Optimization (PSO) and Pareto Dominance [START_REF] Reyes-Sierra | Multi-Objective Particle Swarm Optimizers: A Survey of the State-of-the-Art[END_REF], in order to complete the overall design optimization loop. This algorithm has been implemented in a Matlab environment; the weighting factor was set at 0.7, and the social learning factor and cognitive learning factor both at 1. Leaders were selected from an elitist archive. A crowding factor was used to limit its size to 200.

The leader choice strategy is described in Fig. 6. All of the swarm was evenly spread across the archive so as to assign a leader to each particle. The first particles of the swarm thus always minimize the first objective, while the last particles optimize the second objective. The size of the swarm has been set at 2,000, and 500 generations were chosen. Total computation time lasted roughly 1,000 s on a quad-core PC. To ensure good convergence, the algorithm was run 10 times. The results discussed below represent a combination of the results from each run.

B. Application: Optimization of an all-electric

conversion chain for a wave energy converter SEAREV (French acronym for Autonomous Electrical System of Wave Energy Conversion) is a wave energy conversion system [START_REF] Babarit | Hydrodynamic optimisation and optimal control of a Wave Energy Control[END_REF][START_REF] Ruellan | Sizing methodology of a wave energy converter[END_REF] composed of a float actuated by swell movement. Inside this device, an eccentric mass wheel oscillates relative to the float. An electric generator operating in direct drive converts the mechanical energy of the pendular wheel into electrical energy, thereby causing damping of the relative motion between these two parts. One strategy for controlling the conversion torque consists of generating, via the associated machine drive, a viscous friction torque, e.g. introducing a resistant torque proportional to the angular velocity in addition to a leveling of the converted mechanical power. For each operating point, torque can thus be correlated with angular velocity Ω by [START_REF] Figueroa Barnier | Modélisation des entraînements à grande plage de vitesse en vue de leur conception[END_REF].

Ω Ω Ω Ω (18)
The operating profile considered in this work was obtained for a viscous damping coefficient equal to 5.10 . . and a clipping power of 1 , both having been optimized for the state of an ocean swell 3 , 8 . This temporal angular velocity and torque profiles are depicted in Fig. 8. The time step was set at 0.1 s with a total cycle duration of 1,200 s. The total number of operating points thus equals 12,000. In order to reduce the number of points to be calculated, speed values were discretized and the probability of occurrence for each of the resulting discrete values was then calculated. Fig. 9 shows the simplified cycle placed within the torque-speed plane. The area of each plotted point is proportional to the probability of appearance times the mechanical power ( Ω ∑ ⁄ . The SEAREV system therefore provides an application example of this optimization method on a cycle running in energy generation mode. The level of power being considered allows adopting the hypothesis previously described on transient thermal states, with the temperature rise constraint then being calculated on the basis of the average value of copper and iron losses. Only the converter on the machine side has been included in this example when calculating costs and losses. The scale law coefficients on the power electronic components were derived for a maximum voltage rating of 3,300 V.

In this example, the direct axis current is optimized at each operating point in order to minimize total machine drive losses. The numerical method employed to perform this constrained optimization has been described in the section III. The number of discrete values tested for has been fixed at 100 within the interval ; 0 . The optimization results are presented in Fig. 11 and12. The Pareto front indicates the set of best solutions relative to both the electric chain cost and average electric power , which is a reflection of the produced electrical energy (e.g. the mechanical power (negative) at which the sum of losses (positive) was added). The distribution of losses and costs for the four specific solutions is displayed in Fig. 12. The parameters of these four solutions are then detailed in Table III

VI. CONCLUSION

This paper has presented a machine drive optimization methodology applied on an operating profile with a strong coupling between the machine and the power electronic converter. The aim was optimizing the design of an allelectric conversion chain for a wave energy conversion system. This methodology proves to be applicable to any system in which a machine drive is operating at variable A number of simplifying hypotheses were adopted, primarily to save computation time, yet the integration of more complex models would not pose a problem from a methodological standpoint. A machine model based on finite elements could even be used to more accurately calculate magnetic losses, particularly under conditions of considerable flux weakening, although computation time would rise substantially. A transient thermal model could also improve quality of results especially for applications where thermal time constants of the machine are near from the loss variation periods. VII.

  Fig. 1: Geometric parameterization of the PM machine

Fig. 2 :

 2 Fig. 2: Constrained optimization of , for minimizing

Fig. 4 :

 4 Fig. 4: Equivalent thermal circuit of the PM machine

Fig. 5 :

 5 Fig. 5: Operating cycle design -The optimization flowchart

Fig. 6 :

 6 Fig. 6: Leader choice strategy for the considered MOPSO algorithm

Fig. 7 :

 7 Fig. 7: Schematic diagram of the SEAREV WEC operating principle

Fig. 8 :Fig. 9 :

 89 Fig. 8: Operating profiles (angular velocity and torque) vs. time -Total cycle time equals 1,200 s

  , and Fig 10 shows the operating area capable of being reached by Solution A along with the simplified operating profile. The shading indicates the optimized constrained leading angle (in generator mode Ψ tan , / , ) at each operating point for the purpose of minimizing machine drive losses.

Fig. 10 :Fig. 11 :Fig. 12 :

 101112 Fig. 10: Contour plot of the optimized leading angle in degrees in the (P,Ω) plane for the solution A -Here, the area of each point is proportional to its frequency of appearance

  for electric traction or renewable energy conversion.

  .

	Bi-Objective Sizing Optimization of a PM
	Machine Drive on an Operating Profile
	J. Aubry, H. Ben Ahmed, B. Multon
	S
	(2)

  in Table I below.

TABLE III FOUR

 III DETAILED SOLUTIONS OF THE PARETO FRONT

		A	B	C	D
	P average -230 kW	-234 kW	-238 kW	-240 kW
	C drive	180 k€	250 k€	400 k€	600 k€
	S conv 1.1 MVA 1.3 MVA	2.1 MVA	2.9 MVA
	e	5.1 mm	4.8 mm	5.4 mm	5.9 mm
	h mag	11 mm	15 mm	18 mm	20 mm
	h rot	3 cm	4 cm	7 cm	9 cm
	h slot	18 cm	28 cm	35 cm	42 cm
	h stat	3 cm	6 cm	10 cm	16 cm
	k cf	56%	51%	52%	55%
	p	126	80	51	36
	l u	55 cm	51 cm	65 cm	82 cm
	r s	4.8 m	4.6 m	4.5 m	4.4 m
	n s	378	400	357	342
	M machine	25 t	41 t	70 t	110 t
	L m	3.3 mH	6.4 mH	13 mH	22 mH
	L l	4 mH	7.5 mH	9.2 mH	11 mH
	R s	24mΩ	20mΩ	17mΩ	15mΩ
	Φ e	7Wb	12Wb	21Wb	31Wb

This work was supported by the Bretagne regional council under Grant ARED SEAREV_2.


Mathieu Beau

Valentin A Zagrebnov

The second critical point for the Perfect Bose gas in quasi-one-dimensional traps

Keywords: numbers: 05, 30, Jp, 03, 75, Hh

We present a new model of quasi-one-dimensional trap with some unknown physical predictions about a second transition, including about a change in fractions of condensed coherence lengths due to the existence of a second critical temperature Tm < Tc. If this physical model is acceptable, we want to challenge experimental physicists in this regard.

1. It is well known since 1925 [START_REF] Einstein | Physikalisch-mathematische Klasse[END_REF] for an ideal gas of identical bosons of mass m in a three dimensional cubic box Λ = L 3 described the grand canonical ensemble (V, T, µ) with T the temperature and µ the chemical potential, exists a critical density of thermal gas saturation of ρ c = ζ(3/2)/λ 3 β where λ β = 2πβ/m is the thermal length and β = 1/k B T . The hypothesis of A. Einstein back in 1938 [START_REF] London | [END_REF] by F. London assumed that for a particle density ρ above the critical density ρ c surplux the particles ρ -ρ c starts in the ground state mode k = 0 . Thus we note that ρ k=0 = ρ -ρ c for ρ > ρ c . But this argument macroscopic population of the ground state depends drastically on the geometry of the box. This is only first discovered H. Casimir [START_REF] Casimir | On Bose-Einstein condensation[END_REF] in 1968 for a certain geometry anisotropic and what has motivated the work of M. van den Berg, J. Lewis and J. Pule from 1981 [4] - [7]. These are generalized the concept of Bose-Einstein condensate defining density in thermodynamic limit for an ideal Bose gas in a rectangular box any Λ = L x ×L y ×L z :

ρ 0 = lim η→0 lim Lx,Ly,Lz→∞ k η ρ k = ρ -ρ c . (1)
This definition includes all possibilities from a condensate, given here from [6] and revisited [8] in reference to new terminologies [9]:

(i) The condensate conventional or usual [8], [9] formed by the ground state (type I [6]).

(ii) The fragmented condensate [8], [9] if the generalized condensate is distributed over a set (finite / infinite) modes macroscopically occupied in an energy band near the fundamental mode ie [6]).

N 0 = k kc N k , with N k = O(N ) (Type I / II
(iii) The quasi-condensate [8], [9] if the condensate is distributed over a generalized set of modes occupied mesoscopic in an energy band near the fundamental mode ie N 0 = k kc N k , with N k = O(N δ ), δ < 1 (type III [6]).

In 1983 [5], van den Berg proposed a box model threedimensional anisotropic exponentially in both directions to define a second critical density ρ m separating a regime of condensate generalized type III (for ρ c < ρ < ρ m ) towards a regime of generalized type I condensate (for ρ > ρ m ). Lately, we have reexamined this model [10] questions that we call two-dimensional and we have shown that this a transition between two regimes: quasicondensate and coexistence between a regime of conventional condensate and quasi-condensate, notably by showing that the second critical density ρ m corresponds to a saturation density of gas near condensation similarly to the saturation of thermal gas for ρ c . In addition we have calculated the second critical temperature and determined the fractions similar condensed changed. Then we calculated the effects of this transition on the coherence length for these two regimes of condensate generalized.

However this transition in the model boxes van den Berg is not the case for almost one-dimensional (quasi-1D). Our finding presented here is the existence of this transition between two regimes for the ideal gas of bosons in a harmonic trap quasi-1D. The purpose of our letter is to show where is the transition between a quasi-BEC and BEC conventional find the geometric model for the trap to his observation and to characterize this transition (Fractions condensed, coherence lengths). This theoretical prediction is interesting because we have seen in the experimental [11] and theoretical articles [12], [13] the existence of a second critical temperature called T φ for quasi-1D environments, but calculated and explained differently from ours. 2. Consider an ideal gas of bosons in a harmonic trap characterized by three-dimensional pulsation ω x , ω y , ω z : where external potential is given by:

V ext (r) = 1 2 mω 2 x x 2 + 1 2 mω 2 y y 2 + 1 2 mω 2 z z 2 . ( 2
)
The Hamiltonian for a particle is given by [14]:

T (N =1) Λ = - 2 2m ∆ + V ext (r) . (3)
The energy levels are then given by:

ǫ sxsysz = ω z (s z + 1 2 ) + ω z (s z + 1 2 ) + ω z (s z + 1 2
). ( 4)

Let µ the chemical potential of this system in the grand canonical ensemble, its value must be strictly lower than the ground state to ensure the convergence distribution the Bose-Einstein statistics showing the number average boson in each of the modes of expression N sxsysz = (exp β(ǫ sxsy sz -µ) -1) -1 . Thus we can define the effective chemical potential relative to ground level as µ = ǫ 000 -∆µ. When the potential ∆µ > 0 tends to 0, the thermal gas is saturated and we obtain a critical number of particles when ω 0 = (ω x ω y ω z ) 1/3 tends to 0:

N ≃ 1 ( βω 0 ) 3 R 3,+ d 3 ǫ 1 e β(ǫ+∆µ) -1 < ζ(3) ( βω 0 ) 3 = N c . ( 5
)
The of the spectrum for this trap is different from that of a particle in a box, that is why the type of geometry for the second transition will be different. If we consider that this trap is anisotropic with a quasi-1D type "cigar" in the axis z, we must have the condition ω z ≪ ω ⊥ with ω x = ω y = ω ⊥ . Usually [14], [12], [13] is taken as a model ω ⊥ = λω z with the condition λ ≫ 1. Thus we have an anisotropy level of modes in each direction for ω z ≪ ω ⊥ . If in this case we have an effective chemical potential ω z ≪ ∆µ ≪ ω ⊥ then we can calculate the number of atoms in a condensate generalized (quasi-condensate) in an energy band near the fundamental on the axis z:

N qbec ≡ s=(0,0,sz ) N s ≃ 1 βω z R + dǫ z 1 e β(ǫz+∆µ) -1 ≃ - 1 βω z ln[β∆µ] = N -N c . (6)
Thus for a macroscopic number of particles in the condensate, we must have condition ∆µ = β -1 e -βωz(N -Nc) + . . .. But there is another condition for completing the energetic part along the axis z because the chemical potential ∆µ determines the size of the energy band, decreases when the number of particles N increases and therefore a second critical number of particles noted N m , the condition ω z ≪ ∆µ is false. For this model unlike gas trapped in a box, there is no global notion of "density" of particles as the containment takes place throughout the three-dimensional space R 3 . However we can introduce a scaling for the number of particles n ≡ N ω 3 0 which converges when N tends to infinity and when ω 0 tends to 0 [14]. We then find that when ω ≃ ∆µ is a β ω z ≃ e -β(nm-nc)/ω 2 ⊥ . To resolve this condition, we introduce the model of anisotropic quasi-1D trap following ω z = ω ⊥ e -ω 2 c /ω 2 ⊥ where ω c > 0 is a pulse of anisotropy exponentially. Thus there is a limit to the number of particles in the condensate distributed over the energy axis z:

N m = N c + ω 2 c βω 3 0 . ( 7
)
FIG. 1: Illustration of phase diagram, see ( 5), ( 7)

Now if we go beyond this critical number N m , we shall complete the ground state:

N bec ≡ N s=0 = 1 e β∆µ -1 ≃ 1 β∆µ = N -N m , (8)
and thus the chemical potential is ∆µ = 1/(β(N -N m )) and the value of the quasi-condensate is saturated:

N Qbec ≃ - 1 βω z ln[β∆µ] = N m -N c , since N -N c = (n -n c )/ω 3 0 and thus ω 3 0 N Qbec ≃ ( β) -1 ω 2 ⊥ ln[ω 3 0 /(n -n c )] ≃ ω 2 c ( β) -1
Note that this model of exponential anisotropy is necessary to limit "thermodynamics", ie when ω ⊥ tends to 0, to define the second critical number of particles N m which is the same order as the first critical number N c (proportional to ω -3 0 ). For finite pulses in a given experimental situation, we must calculate the parameter ω c depending on the aspect ratio λ = e ω 2 c /ω 2 ⊥ and see if the ratio ω c /ω ⊥ is sufficiently large to consider that one is in an exponential phase. Usually if one considers a system where very anisotropic λ ≫ 1 to consider this point of view that the system is exponentially anisotropic, we have to have ln(λ) ≫ 1 .

3. In experiments with BEC, it is important to know the critical temperatures associated with corresponding critical densities. The first critical temperatures: T c , for a given particles number N is given by solutions of this equation: N = N c = ζ(3)/( βω 0 ) 3 . Then we get:

T c = ω 0 k B ζ(3) 1/3 N 1/3 . ( 9
)
By the expression for the second critical particles number one gets the following relation between the first and the second critical temperatures because T m are given by solution of this equation N = N m = N c + ω 2 c /( βω 3 0 ): Here τ = ω c /k B ζ(3) 1/2 is "effective" temperatures related to the corresponding geometrical shape. We can compute the relative difference between both critical temperature due to (10) , see figure 2:

T 3 m + τ 2 T m = T 3 c . ( 10
)
T c -T m T c = f ( N ωc N ) (11)
where

N ωc = ω 3 c ζ(3) 1/2 ω 3 0 defined by Nω c N = τ 3 T 3 c
and with:

f (x) = 1 -Ω(x) 1/3 + x 2/3 3 Ω(x) -1/3 ( 12
)
where Ω(x) = 1 2 (1 + 1 + 4 27 x 2 ). The temperature dependence of the "quasicondensate" N Qbec is, see figure 3:

N Qbec N =    1 -T Tc 3 , T m ≤ T ≤ T c , τ 2 T 3 c T , T ≤ T m . ( 13
)
The corresponding ground state conventional BEC is, see figure 3:

N bec N = 0 , T m ≤ T ≤ T c , 1 -( T Tc ) 3 1 + τ 2 T 2 , T ≤ T m , (14)
and for the two coexisting condensates one gets:

N -N c N = N 0 N = N Qbec + N bec N = 1 - T T c 3/2 . (15)
4. Another physical observable to characterize this second critical temperature is the correlation function condensate coherence length. In quasi-1D harmonic trap, the correlation function is locally defined by:

g(r, r ′ ) = s φ s (r)φ s (r ′ ) e β(ǫs-µ) -1 , (16)
where the wave function φ s (r) are the well known three dimensional hermite functions:

φ s (r) = φ sx (x)φ sy (y)φ sz (z),
with for ν = x, y, z (r x ≡ x, r y ≡ y, r z ≡ z):

φ sν (r ν ) = c sν 1 √ L ν H sν (r ν /L ν )e - r 2 ν 2L 2 ν , with c sν = 1/ s ν !2 sν √
π, ν = x, y, z and where the effective confinement lengths L ν , ν = x, y, z are given by: L ν = mων . Notice that the limiting diagonal function n(r) ≡ σ(r, r) is local space particle density at point r.

To compute the correlations in the direction z, induced by the quasi-condensate, we have to consider r = r -∆r/2 r ′ = r + ∆r/2, where ∆r/2 = (0, 0, ∆z) and integrate r over the space to keep the effect of condensate through the trap, it give the mean value of correlation for the distance ∆r: = g N bec (∆r) + g Qbec (∆r) + g bec (∆r) [START_REF] Beau | [END_REF] where correlation function is divided in three parts, each part is the truncation of the sum ( 16) : the noncondensate correlations g N bec ≡ sx>0,sy>0,sz , the quasi-condensate correlations g Qbec ≡ sx=0,sy=0,sz and the condensate correlations g bec ≡ (correlation of the ground state). By formulas [START_REF] Abramowitz | Handbook of Mathematical functions with formulas, graphs and mathematical tables[END_REF] and [START_REF] Beau | [END_REF] and by mathematical assumptions [START_REF] Abramowitz | Handbook of Mathematical functions with formulas, graphs and mathematical tables[END_REF], [START_REF] Beau | [END_REF] for the noncondensate part, we get, see figure 4:

g (∆r) =
g N bec (∆r) ≃ R 1 d x R 1 d y R 1 d z ∞ j=1 e jβµ πj ω 0 /2 3 × e -jmβ(ω 2 x x 2 +ω 2 y y 2 +ω 2 z z) 2 e -2m (∆r) 2 jβ 2 = ∞ j=1 e jβµ 1 j 3 3 ω 3 0 exp[-2m (∆r) 2 jβ 2 ], (18)
which decrease exponentially fast for µ < 0. The quasicondensate part, for T < T c , by the scaling argument of localization of the scale of energetic momentum of quasicondensate given by ∆µ, and by mathematical assump- tions [START_REF] Abramowitz | Handbook of Mathematical functions with formulas, graphs and mathematical tables[END_REF], [START_REF] Beau | [END_REF], is given by, see figure 4:

g Qbec (∆r) ≃ R 1 d z ∞ j=1 e jβ∆µ 1 πj ω 0 /2 × exp[-jmβ(ω 2 z z) 2 β∆µ] exp[-2m (∆r) 2 β∆µ jβ 2 ], = ∞ j=1 e jβµ 1 j ω z exp[-2m (∆z) 2 ∆µ jβ 2 ]. (19)
The condensate part, for T < T m is given by:

g bec (∆r) = R 3 d rφ 000 (∆r)φ 000 (-∆r)N bec = N bec exp ( ∆z 2 L 2 z ) ≃ N bec . (20)
We can defined the coherence length of the quasicondensate as the maximal scale of length of ∆z such that the correlations are non-negligible when ω z , ω ⊥ tends to zero. Then by a scaling argument, see equation (19) we obtain L c = /mω ⊥ β∆µ, where ∆µ is given by β -1 e βωz(N -Nc) for N c < N < N m . Thus L c mω ⊥ / = (ω ⊥ /ω z ) γ(T ) , where the exponent γ(T ) = βω . This gives us a relative difference between the two critical densities of 0.587, see relation 11 and figure 2 with the implementation of the interactions between particles and its effect on the second critical temperature as it is done for the first critical temperature [15]. An other question for quasi-two-dimensional boxes is the validity of Hohenberg theorem [9] which forbidden the existence of the condensate on the ground state for homogeneous system.

FIG. 2 :FIG. 3 :

 23 FIG. 2: Curve illustrating (11) as a function of τ Tc .

R 3 d 3

 33 rg( r -∆r/2, r + ∆r/2).

15 FIG. 4 :

 154 FIG.4:(a) Allure of correlation functions (19) with ω0 = 1, β 2 /2m = 1 and β∆µ = 1, 0.1, 0.01, 0.001, 0.00001, 0.0000001 respectively curves (---), (---), (...) and (-). Notice that (...) and (-) are combined, what can be interpreted as a saturation of the coherence of the thermal gas. (b) Allure of correlation functions (18) with ω0 = 1, β 2 /2m = 1 for β∆µ = 0.001, 0.0001, 0.000001 respectively (---), (---), (-).

3 - 1 ,) 7 .FIG. 5 :

 3175 FIG. 5: curve illustrating the exponent (21) as a function of T Tc and with Tc τ = 2 3. This gives us a relative difference between the two critical densities of 0.587, see relation 11 and figure2

We are thankful to Markus Holtzmann for useful discussions about the exponential anisotropy effect.


Guillaume Merle

Jean-Marc Roussel

Jean-Jacques Lesage

Improving the Efficiency of Dynamic Fault Tree Analysis by Considering Gates FDEP as Static

This paper focuses on one of the dynamic gates which are used in Dynamic Fault Trees (DFT), which is the Functional Dependency (FDEP) gate. Gate FDEP has been considered as equivalent to a set of OR gates in the literature, but this equivalence has seldom been exploited for the analysis of DFTs. In this paper, we show that in most cases, starting from a DFT including FDEP gates, the use of this static equivalence provides significant advantages for DFT Analysis.

INTRODUCTION

Fault Tree Analysis (FTA) is one of the oldest, most diffused techniques in industrial applications, for the dependability analysis of large safety-critical systems [START_REF] Henley | Reliability Engineering and Risk Assessment[END_REF][START_REF] Leveson | Safeware: System Safety and Computers[END_REF][START_REF] Stamatelatos | Fault Tree Handbook with Aerospace Applications[END_REF]. FTA is usually carried out at two levels: a qualitative level in which the list of all the possible combinations of events that lead to the Top Event (TE) is determined (the minimal cut sets); and a quantitative level in which the probability of the occurrence of the TE, and of the other nodes of the tree, is calculated. One of the main restrictive assumptions in FTA is that basic events must be assumed to be statistically independent, and their interaction is described by means of boolean OR/AND gates, so that only the combination of events is relevant, and not their sequence. We refer to this model as Static Fault Tree (SFT). Dugan et al. [START_REF] Dugan | Dynamic fault-tree models for fault-tolerant computer systems[END_REF][START_REF] Dugan | Developing a low-cost high-quality software tool for Dynamic fault-tree analysis[END_REF] proposed a new model allowing to include various kinds of temporal and statistical dependencies in the SFT model, which is the Dynamic Fault Tree (DFT). The DFT is based on the definition of new gates: Priority-AND (PAND), Functional Dependency (FDEP), Warm Spare (WSP), and Sequence Enforcing (SEQ). Gate FDEP was introduced in [START_REF] Dugan | Dynamic fault-tree models for fault-tolerant computer systems[END_REF] and allows to model common cause failures by using the concept of preemption: a trigger event allows to force dependent basic events to fail (independently of their own failure). Even though gate FDEP is most often considered as a dynamic gate [START_REF] Boudali | Dynamic Fault Tree analysis through input/output interactive Markov chains[END_REF][START_REF] Stamatelatos | Fault Tree Handbook with Aerospace Applications[END_REF], some authors have already mentioned that it is equivalent to a set of OR gates [START_REF] Ejlali | FPGA-based Monte Carlo simulation for fault tree analysis[END_REF]. However, to the best of our knowledge, the influence of this equivalence on DFT analysis has never been studied. In the present paper, thanks to a DFT example, we show the influence of this equivalence on DFT Analysis. This study has been developed under the following assumptions:

• we focus on modular approaches, for which the static and dynamic parts of a DFT are analyzed by means of specific dedicated approaches;

• the results presented hold for approaches based on Continuous Time Markov Chains as well as Stochastic Petri Nets;

• this study has been validated by the tool Galileo [START_REF] Dugan | Developing a low-cost high-quality software tool for Dynamic fault-tree analysis[END_REF].

The definition of gate FDEP is recalled in Section 2. The advantages of the equivalence between gate FDEP and a set of OR gates for DFT Analysis are studied in Sections 3 to 5.

DEFINITION OF GATE FDEP

According to [START_REF] Boudali | Dynamic Fault Tree analysis through input/output interactive Markov chains[END_REF][START_REF] Dugan | Dynamic fault-tree models for fault-tolerant computer systems[END_REF][START_REF] Stamatelatos | Fault Tree Handbook with Aerospace Applications[END_REF], the FDEP gate -Functional Dependency gate -is a dynamic gate comprised of a trigger input event -either a basic event or the output of another gate of the tree -and a set of dependent basic events. Fig. 1 provides a pictorial depiction of an FDEP gate with 2 dependent basic events A and B, T representing the trigger event. When the trigger event occurs, the dependent basic events are forced to Numerous examples of the use of gate FDEP can be found in the literature which show that, in particular, the semantics of this gate is well adapted to model common cause failures. Indeed, the FDEP gate in Fig. 1 can allow to model a system with a component T whose failure is a common cause failure leading to the simultaneous failures of two basic events A and B, keeping in mind that basic events A and B can also fail by themselves. Gate FDEP thus offers the engineers an easy way to translate such preemptive behaviors in DFTs.

Some authors have suggested that gate FDEP could be considered as a set of OR gates. For example, the authors of [START_REF] Stamatelatos | Fault Tree Handbook with Aerospace Applications[END_REF], who proposed gate FDEP as a dynamic gate, also provide an alternative static representation which only includes OR gates -one for each basic event. Besides, in [START_REF] Ejlali | FPGA-based Monte Carlo simulation for fault tree analysis[END_REF], the time to failure of the dependent components of the FDEP gate is defined as equal to the minimum of the time left before the occurrence of the trigger event and the time left before the occurrence of the dependent event. This corresponds to the temporal definition of operator OR, based on the dates of occurrence of its input events, as described in [START_REF] Merle | Algebraic Expression of the Structure Function of a subclass of Dynamic Fault Trees[END_REF][START_REF] Merle | Probabilistic Algebraic Analysis of Fault Trees With Priority Dynamic Gates and Repeated Events[END_REF].

3 MODULAR ANALYSIS OF DYNAMIC FAULT TREES We are going to illustrate the usefulness of the equivalence between gate FDEP and a set of OR gates for the qualitative and quantitative analysis of DFTs thanks to the example of a computer system from (Stamatelatos and Vesely 2002) shown in Fig. 3. This DFT models the potential failure of the system in Fig. 2.

Figure 2: System corresponding to the DFT in Fig. 3 The system in Fig. 2 is composed of:

• a redundant bus subsystem including two identical buses, of which one is required for system operation;

• a redundant processing subsystem including two redundant processors A1 and A2 and a spare processor A, which can replace either upon failure.

The subsystem can continue to operate until all three processors have failed;

• a redundant memory subsystem including five memory units -M 1 to M 5 -of which three are required, these memory units being connected to the redundant bus via two memory interface units, memory unit 3 -M 3 -being connected to both interfaces for redundancy;

• an application subsystem including a human operator, an application (software (SW)), and a graphical user interface (hardware (HW)).

A complete description of the architecture, as well as the explanation of each considered failure, can be found in [START_REF] Stamatelatos | Fault Tree Handbook with Aerospace Applications[END_REF].

The DFT shown in Fig. 3 can be decomposed into 4 subtrees [START_REF] Dugan | Developing a low-cost high-quality software tool for Dynamic fault-tree analysis[END_REF] according to the gates that it contains. Subtrees 3 and 4 contain static gates, only, so they can be considered as static. However, subtree 1 contains Spare gates and is hence dynamic. Finally, subtree 2 contains gates FDEP, 3-out-of-5, and AND, so it can be considered as static or dynamic, depending on the model considered for gate FDEP.

Many researchers have explored the use of divideand-conquer approaches for analyzing such DFTs [START_REF] Bobbio | Parametric Fault Trees with Dynamic Gates and Repair Boxes[END_REF][START_REF] Chatterjee | Modularization of fault trees: A method to reduce cost of analysis[END_REF][START_REF] Dugan | Developing a low-cost high-quality software tool for Dynamic fault-tree analysis[END_REF][START_REF] Dutuit | A linear-time algorithm to find modules in fault trees[END_REF][START_REF] Rosenthal | Decomposition methods for fault tree analysis[END_REF]). For instance, in [START_REF] Dugan | Developing a low-cost high-quality software tool for Dynamic fault-tree analysis[END_REF], a subtree is marked as dynamic, and solved by using Markov Chains, if it contains at least one dynamic gate. If a subtree contains no dynamic gates, it is classified as static and solved using BDD-based methods. The tool Galileo [START_REF] Dugan | Developing a low-cost high-quality software tool for Dynamic fault-tree analysis[END_REF] is based on this method. The authors of [START_REF] Bobbio | Parametric Fault Trees with Dynamic Gates and Repair Boxes[END_REF]) also exploit the concept of modularity, a subtree that is statistically independent from the rest of the FT being denoted as a module: static modules can be analyzed by means of suitable combinatorial techniques whereas dynamic modules require a state-space analysis which is obtained by translating the dynamic module into a Generalized Stochastic Petri Net.

In both cases (Markov Chains vs. Stochastic Petri Nets), it can be useful to be able to reduce the dynamic Figure 3: The 4 subtrees of the DFT in [START_REF] Stamatelatos | Fault Tree Handbook with Aerospace Applications[END_REF] part of the DFT because the computational complexity of its solving is considerably higher than the computational complexity of the solving of the static part of the DFT.

The impact of the use of the equivalence between gate FDEP and a set of OR gates on the qualitative and quantitative analysis of this DFT example is analyzed in Sections 4 and 5, respectively.

4 ADVANTAGES FOR QUALITATIVE ANALY-SIS When an independent subtree of a DFT is static, BDD-based methods allow to determine the minimal cut sets of the DFT [START_REF] Coudert | Fault Tree Analysis: 10 20 prime implicants and beyond[END_REF][START_REF] Dutuit | Exact and truncated computations of prime implicants of coherent and noncoherent fault trees with Aralia[END_REF]. However, when such an independent subtree of a DFT is dynamic, the concept of minimal cut set must be extended to the concept of minimal cut sequence [START_REF] Tang | Minimal Cut Set/Sequence Generation for Dynamic Fault Trees[END_REF], and BDDs can no longer be used to perform the qualitative analysis of such FTs. The authors of [START_REF] Tang | Minimal Cut Set/Sequence Generation for Dynamic Fault Trees[END_REF]) address this problem by using specific BDDs, denoted as Zero-suppressed BDDs, to perform the qualitative analysis of both the static and dynamic independent subtrees of DFTs. Other authors, such as [START_REF] Bobbio | Parametric Fault Trees with Dynamic Gates and Repair Boxes[END_REF], convert the dynamic subtree into a Generalized Stochastic Petri Net whose occurrence graph allows to determine the minimal cut sequences of the DFT.

The fact to consider gate FDEP as dynamic or static does not have much impact on the methods used to perform the qualitative analysis of DFTs, but on the solution obtained. Indeed, in some cases, a set of minimal cut sequences can be equivalent to a single minimal cut set: for instance, if A and B are two non-repairable basic events, the set of two minimal cut sequences {[A, B] , [B, A]}1 is equivalent to the minimal cut set A • B. As a consequence, if gate FDEP is considered as dynamic, it will provide minimal cut sequences which will be equivalent to the less numerous minimal cut sets that it will provide if it is considered as static.

Let us illustrate this possible simplification on the DFT in Fig. 3. If we consider gate FDEP as equivalent to a set of OR gates, subtree 2 contains 1 3-out-of-5 gate, 1 AND gate and 3 FDEP gates, and can hence be considered as static and equivalent to the SFT in Fig. 4.

Figure 4: SFT equivalent to the subtree 2 of the DFT in Fig. 3 The subtree 2 in Fig. 3 has 74 minimal cut sequences whereas it has only 17 minimal cut sets in Fig. 4. Both results are equivalent, since a single cut set may represent a more or less big set of cut sequences. For instance, the 6 order-3 min-

imal cut sequences [M 1, M 2, M 3], [M 1, M 3, M 2], [M 2, M 1, M 3], [M 2, M 3, M 1], [M 3, M 1, M 2], and [M 3, M 2, M 1] are equivalent to the minimal cut set M 1 • M 2 • M 3.
Even though both results are equivalent, minimal cut sets represent a more concise -and hence more useful -result to the practitioner than the corresponding set of minimal cut sequences. It can hence be interesting for the qualitative analysis of DFTs to consider gate FDEP as equivalent to a set of OR gates.

This equivalence also provides some advantages for the quantitative analysis of DFTs, as explained in Section 5.

ADVANTAGES FOR QUANTITATIVE ANAL-

YSIS The quantitative analysis of DFTs underlies the numerical computation of stochastic models, and consequently the use of specific software tools. For a given software, evaluation criteria are needed to compare the computational complexity for the determination of the failure probability of the TE for both models of gate FDEP. These criteria are studied in the following sections.

Experimental protocol

We chose to use the well-known DFT analysis tool Galileo [START_REF] Dugan | Developing a low-cost high-quality software tool for Dynamic fault-tree analysis[END_REF], which exploits the concept of modularity (BDDs for static subtrees and Markov Chains for dynamic subtrees).

Since Galileo, as most DFT analysis softwares, is dedicated to personal computers, its performances are hence limited by the processor performance and the size of the RAM of the computer, since the algorithms used have an exponential complexity and the constructed state spaces of the models used are stored in the RAM, thus leading to a risk of saturation. We need evaluation criteria to evaluate the performances of the tool Galileo: the evaluation criteria that we chose for this study hence are the computation time and the memory usage. It can be noted that the computation time and memory usage needed to perform the quantitative analysis of subtrees 1, 3, and 4 of the DFT in Fig. 3 will not be impacted by the equivalence between gate FDEP and a set of OR gates since subtree 2 is statistically independent from the other subtrees. As a consequence, the impact of this equivalence on the quantitative analysis of the DFT in Fig. 3 will be limited to the quantitative analysis of subtree 2, since it is the only subtree containing FDEP gates.

Besides, the performances evaluated can be im-pacted by some parameters of the study. As we want to study the impact of the equivalence between gate FDEP and a set of OR gates on the quantitative analysis of DFTs, the first influent parameter coming to mind is the number of FDEP gates, whose impact will be studied in Section 5.3. Besides, it can be interesting to study what the impact of time granularity -the number of failure probabilities calculated from time zero to the mission time -is on quantitative analysis in both cases. The influence of this parameter is studied in Section 5.2.

Impact of time granularity

The computation time and memory usage obtained with Galileo to perform the quantitative analysis of the subtree 2 of the DFT in Fig. 3 when gate FDEP is considered as dynamic or static are presented in Table 1. The quantitative analysis was performed on a Pentium 4 processor with 4 GB RAM for a mission time of T = 10, 000 hours, and with a time granularity varying from 10 to 10, 000. The memory usage of Galileo -approximately 7 MB -was removed from the memory usage obtained to get the memory usage needed for the quantitative analysis in itself.

It can be seen in Table 1 that the fact to consider gate FDEP as equivalent to a set of OR gates allows to divide the computation time for the failure probability of the TE by a factor of between 6 and 15. Even though the memory usage is a linear function of the time granularity in both cases, it can be divided by a factor of approximately 2 by exploiting this equivalence, thus allowing to reach a higher time granularity on the same computer. 3. Subtree 2 models the behavior of the subsystem in Fig. 5. The analysis of this subsystem allows to understand the corresponding DFT: the failure of the Memory In- This subsystem can be generalized to see the impact of the number of FDEP gates on the complexity of quantitative analysis. Let us consider the subsystem in Fig. 5 as the order-2 subsystem since it has 2 Memory Interface Units; its DFT contains 3 FDEP gates. The order-n subsystem with n Memory Interface Units is shown in Fig. 8. In the same way that subtree 2 in the DFT in Fig. 3 describes the potential failure of the subsystem in Fig. 5, the subsystem in Fig. 8 can be described by a DFT which is shown in Fig. 6 and contains (2n -1) FDEP gates. The equivalent SFT is shown in Fig. 7 and contains (2n + 1) OR gates.

The computation time and memory usage needed by Galileo to perform the quantitative analysis of the DFT in Fig. 6 and of the SFT in Fig. 7 are presented in Table 2.

The results obtained in Table 2 show that the computation time varies with the mission time while the memory usage does not. However, once again, the fact to considerate gate FDEP as equivalent to a set of OR gates allows to divide the computation time for the failure probability of the TE by a factor of between 3 and 25, even though this equivalence becomes fruitful for a high number of FDEP gates, only. The results obtained show that this equivalence has a significant influence on the computational complexity for the determination of the failure probability of the TE of a DFT when time granularity increases. However, this influence is somehow limited when the number of FDEP gates increases, and even least than expected, since it becomes fruitful for a high number of FDEP gates, only.

Even if this aspect is less quantifiable, we reckon that this equivalence is more fruitful on the qualitative analysis of DFTs. Indeed, as shown for the subtree 2 of the DFT example in Fig. 3, a large set of minimal cut sequences can happen to be quite difficult to exploit whereas an equivalence with a small set of minimal cuts may be more useful to the practitioner.

In any case, we have shown in this study that the fact to consider gate FDEP as a static gate can allow to convert dynamic parts of a DFT into static parts, thus simplifying both the qualitative and quantitative analysis of the DFT. It can be noted that it is in the case of dynamic subtrees containing gates FDEP only, and no other dynamic gate, that this equivalence has the most significant impact, since it allows to convert whole dynamic subtrees into static -and hence more readily analyzable -subtrees.

Figure 1 :

 1 Figure 1: An FDEP gate with 2 dependent basic events A and B occur.Numerous examples of the use of gate FDEP can be found in the literature which show that, in particular, the semantics of this gate is well adapted to model common cause failures. Indeed, the FDEP gate in Fig.1can allow to model a system with a component T whose failure is a common cause failure leading to the simultaneous failures of two basic events A and B, keeping in mind that basic events A and B can also fail by themselves. Gate FDEP thus offers the engineers an easy way to translate such preemptive behaviors in DFTs.Some authors have suggested that gate FDEP could be considered as a set of OR gates. For example, the authors of[START_REF] Stamatelatos | Fault Tree Handbook with Aerospace Applications[END_REF], who proposed gate FDEP as a dynamic gate, also provide an alternative static representation which only includes OR gates -one for each basic event. Besides, in[START_REF] Ejlali | FPGA-based Monte Carlo simulation for fault tree analysis[END_REF], the time to failure of the dependent components of the FDEP gate is defined as equal to the minimum of the time left before the occurrence of the trigger event and the time left before the occurrence of the dependent event. This corresponds to the temporal definition of operator OR, based on the dates of occurrence of its input events, as described in[START_REF] Merle | Algebraic Expression of the Structure Function of a subclass of Dynamic Fault Trees[END_REF][START_REF] Merle | Probabilistic Algebraic Analysis of Fault Trees With Priority Dynamic Gates and Repeated Events[END_REF].

Figure 6 :

 6 Figure 6: DFT of the order-n subsystem

Figure 7 :

 7 Figure 7: SFT of the order-n subsystem

Table 1 :

 1 Computation Time CT (s) and Memory Usage MU (MB) obtained with Galileo when gate FDEP is considered as dynamic and equivalent to a set of OR gates Impact of the number of FDEP gates Let us consider the DFT in Fig.

	Time granularity	Galileo dynamic CT MU CT MU static
	10	1	3	0	2
	50	4	11	0	6
	100	9	21	0	10
	500	46	99	3	47
	1,000	93	197	7	94
	5,000	577	980	73 465
	10,000	1,293 1,959 220 930
	5.3

Table 2 :

 2 Variation of the Computation Time CT (s) and Memory Usage MU (MB) obtained with Galileo with the number of FDEP gates when gate FDEP is considered as dynamic or static. Two mission times are considered: T = 100 hours and T = 10, 000 hours. Pr{T E} (10, 000) Pr{T E} (100) Pr{T E} (10, 000)

	Order n	Number of FDEP Pr{T E} (100) gates dynamic CT MU CT	MU	CT	MU	static CT	MU
	2	3	0	0.5	0	0.5	0	0.4	0	0.4
	3	5	0	1	0	1	0	0.5	0	0.5
	4	7	0	3	2	3	0	1.5	0	1.5
	5	9	2	8	5	8	0	3	0	3
	6	11	4	13	13	13	0	7	0	7
	7	13	11	24	41	24	0	15	1	15
	8	15	28	41	104	41	2	30	4	30
	9	17	51	66	193	66	11	59	17	59
	10	19	82	100	300	100	25	105	35	105
	11	21	125	148	445	148	35	186	56	186

The notation [A, B] denotes the sequence of failures in which B fails after A has failed.