Ulisses Alves

Maciel Neto

Paulo Sergio Pereira Da Silva

Karine Beauchard
email: karine.beauchard@ens-rennes.fr

Pierre Rouchon
email: pierre.rouchon@mines-paristech.fr

U A Maciel Neto

† P S Pereira Da Silva

§ P Rouchon

H 1 -Control of an ensemble of half-spin systems replacing Rabi pulses by adiabatic following

INTRODUCTION

Most controllability results available for infinite dimensional bilinear systems are related to systems with discrete spectra (see for instance, Beauchard and Coron, [START_REF] Beauchard | Controllability of a quantum particle in a moving potential well[END_REF] for exact controllability results and Beauchard and Nersesyan, [START_REF] Beauchard | Semi-gobal weak stabilization of bilinear schrdinger equations[END_REF] Chambrion, Mason, Sigalotti, and Boscain, [START_REF] Chambrion | Controllability of the discrete-spectrum schrödinger equation driven by an external field[END_REF], Nersesyan [START_REF] Nersesyan | Growth of sobolev norms and controllability of schrödinger equations[END_REF], for approximate controllability results). As far as we know, very few controllability studies consider systems admitting a continuous part in their spectra. In Mirrahimi [START_REF] Mirrahimi | Lyapunov control of a quantum particle in a decaying potential[END_REF] an approximate controllability result is given for a system with mixed discrete/continuous spectrum: the Schrödinger partial differential equation of a quantum particle in an N-dimensional decaying potential is shown to be approximately controllable (in infinite time) to the ground bounded state when the initial state is a linear superposition of bounded states. Ensemble controllability was studied for instance Li and Khaneja [START_REF] Li | Control of inhomogeneous quantum ensemble[END_REF], [START_REF] Li | Control of inhomogeneous quantum ensemble[END_REF] in the context of quantum systems that are described by Bloch equations depending continuously on a finite number of scalar parameters, and with a finite number of control inputs. The goal of ensemble controllability is to simultaneously steer a continuum of systems between states of interest with the same control inputs. Such continuous family of ordinary differential systems sharing the same control inputs can be seen as an interesting example of infinite dimensional systems with purely continuous spectra. The role of Lie algebras in the characterization of ensemble controllability was the main contribution of Li and Khaneja [START_REF] Li | Control of inhomogeneous quantum ensemble[END_REF]. In Beauchard, Coron, and Rouchon [START_REF] Beauchard | Controllability issues for continuous spectrum systems and ensemble controllability of bloch equations[END_REF], these aspects are studied under a functional analysis setting, developed for infinite dimensional systems governed by partial differential equations (see, e.g., Coron [START_REF] Coron | Control and Nonlinearity[END_REF], for samples of these methods). Many results for this ensemble of Bloch equations, that include the discrimination between approximate and exact controllability are obtained in that paper. For instance, it is shown that a priori bounded L 2 -controls are not sufficient to achieve exact controllability, but unbounded controls (containing a sum of Dirac masses) are able to recover controllability. In Beauchard et al. [START_REF] Beauchard | Stabilization for an ensemble of half-spin systems[END_REF] it is shown that the ensemble of Bloch equations is approximately controllable to the south pole of the Bloch Sphere (in the Sobolev space H1 ) in finite time. with unbounded controls.

In this paper one considers the ensemble M(t, ω) ∈ S 2 ⊂ R 3 of Bloch equations:

Ṁ(t, ω) = S(u(t)e 1 + v(t)e 2 + ωe 3 )M(t, ω), (1)
where -∞ < ω * < ω * < +∞, ω ∈ (ω * , ω * ), and {e 1 , e 2 , e 3 } is the canonical basis of R 3 , and S(•) is the (unique) antisymmetric matrix that defines the wedge product, that is,

S(c)v = c ∧ v for all c, v ∈ R 3 . Note that, for c = c 1 e 1 + c 2 e 2 + c 3 e 3 .
one has:

S(c) =   0 -c 3 c 2 c 3 0 -c 1 -c 2 c 1 0   (2)
It must be stressed that u(t) and v(t) are common controls for all the members of the ensemble, and they cannot depend on ω. Such control problems have been addressed since a long time in the NMR community (see [START_REF] Tycko | Broadband population inversion[END_REF], [START_REF] Kobzar | Pattern pulses: design of arbitrary excitation profiles as a function of pulse amplitude and offset[END_REF]). For simplicity the derivative of M with respect to time is denoted by Ṁ, and the derivative of M with respect to ω is denoted by M . To avoid confusion with derivation with respect to ω, the transpose of a matrix A will be denoted by A . We shall consider only the Larmor dispersion (represented by the parameter ω). One does not consider rf-inhomogeneity in this paper. We will state the control problem:

Definition 1 (Control Problem) Given an initial state M 0 ∈ H 1 ((ω * , ω * ), S 2
)), and given ε > 0, choose T f > 0 and construct (bounded) controls u :

[0, T f ] → R and v : [0, T f ] → R in a way that 1 M(T f , •) + e 3 L ∞ ≤ ε.
Beauchard et al [START_REF] Beauchard | Stabilization for an ensemble of half-spin systems[END_REF] has investigated feedback stabilization of the ensemble [START_REF] Beauchard | Controllability of a quantum particle in a moving potential well[END_REF], which possesses a continuous spectra. As in [START_REF] Mirrahimi | Lyapunov control of a quantum particle in a decaying potential[END_REF], the feedback design is based on a Lyapunov functional that is closely related to the norm of the state space (H 1 ), a Banach space. This feedback assures local L ∞ -stabilization of the initial state M 0 (ω) to -e 3 . The control law that was considered in that paper is a sum of a bounded control (u 1 (t), u 2 (t)) with a comb of T -periodic π-Rabi pulses.

The main contribution of the present work is to show that this comb of Rabi pulses may be replaced by an adiabatic following technique, which relies on a priori bounded T -periodic pair of inputs (u(t), v(t)).

The potential practical interest of such stabilization techniques consists in a simple algorithm providing an open loop control steering from the initial state to the final state. This algorithm just consists in the numerical integration of the closed loop system where the control values are recorded at each sample integration time.

The main results are presented in section 2. Computer simulations are presented in section 3. Those simulations have indicated a faster convergence than the original feedback law of [START_REF] Beauchard | Stabilization for an ensemble of half-spin systems[END_REF] , [START_REF] Beauchard | Stabilization of an arbitrary profile for an ensemble of half-spin systems[END_REF]. The shape of the adiabatic pulses are presented in figure 6 (see appendix A). Some proof sketches are presented in the Appendix B.

MAIN RESULTS AND HEURISTICS

The main ingredients of our control strategy are the adiabatic following and the H 1 -stabilizing strategy. One shall begin describing the adiabatic following and its (discontinuous) propagator.

The adiabatic propagator

Consider the adiabatic propagator equation:

Ȧ(t, ω) = S( ū(t)e 1 + v(t)e 2 + ωe 3 )A(t, ω) (3)
where A(t, ω) ∈ SO(3), and A(0, ω) = I, ∀ω ∈ [ω * , ω * ] and the pair ( ū(t), v(t)) is the adiabatic control.

Consider that one applies the adiabatic control

ū(t) = B 1 (t) sin φ (t), (4a) v
(t) = B 1 (t) cos φ (t), (4b)
where φ (t) and B 1 (t) are defined by:

φ (t) = k(t) ā(t), φ (0) = 0 B 1 (t) = k(t) b(t)
where ā(•), b(•), and k(•) are T -periodic functions defined by ā(t) = a(t/T ), b(t) = b(t/T ), and k(t) = Kk(t/T ), where a(•), b(•), and k(•) are 1-periodic normalized functions defined in the Appendix A, and K > 0 is a chosen gain. We shall re-initialize the propagator to the identity at t 0 k = kT , for k ∈ N. One will denote the discontinuous 2

but T -periodic propagator

A : R × [ω * , ω * ) → SO(3) by A(t, ω). The left-limit lim t→T - A(t, ω) is denoted by A(T -, ω)
The following convergence result can be proved by using standard adiabatic techniques [START_REF] Teufel | Adiabatic Perturbation Theory in Quantum Dynamics[END_REF] or by averaging methods [START_REF] Maciel Neto | Controle adiabático de ensembles quânticos via método das médias[END_REF]:

Theorem 1 Fix K > max{|ω * |, |ω * |}.
Then lim

T →∞ A(T -, ω) -I L ∞ = lim T →∞ max ω∈[ω * ,ω * ] A(T -, ω) - I = 0.
It must be stressed that the last theorem only holds because of the particular symmetries of the adiabatic pulses described in Appendix A (see Figure 6).

The H 1 control law of the auxiliary system

One forgets for the moment that A(t, ω) has some discontinuities and considers the following rotating coordinate change:

N(t, ω) = A (t, ω)M(t, ω).
In each interval [kT, (k + 1)T ) one may show that 3 :

Ṅ(t, ω) = S A (t, ω)(u 1 (t)e 1 + u 2 (t)e 2 ) N(t, ω) (5)
where

u(t) = u 1 (t) + u(t) and v(t) = u 2 (t) + v(t) and N(kT, ω) = A (kT, ω)M(t, ω) = M(t, ω), k ∈ N.
This construction would imply that N(t, ω) may have discontinuities at t = kT, k ∈ N. In this paper one regards the continuous solution N(t, ω) of ( 5):

Remark 1 (Very important) The system [START_REF] Beauchard | Semi-gobal weak stabilization of bilinear schrdinger equations[END_REF] with input (u 1 (t), u 2 (t)) and initial condition N 0 = M 0 will be called by Auxiliary System. Its (continuous) solution will be denoted by N(t, ω). This implies that the discontinuous map M 1 (t, ω) = A (t, ω)N(t, ω) will not be the corresponding (continuous) solution of (1), with input u(t) = u 1 (t) + u(t) and v(t) = u 2 (t) + v(t) and initial condition M 0 . In appendix B it is shown that M 1 (t, ω) = A (t, ω)N(t, ω) is only an approximation of the (continuous) solution M(t, ω).

Consider the Lyapunov functional

L = 1 2 N +e 3 2 H 1 = ω * ω * 1 2 N , N + 1 + N, e 3 dω
In order to compute L note that ξ = u 1 (t)e 1 + u 2 (t)e 2 does not depend on ω. One has

Ṅ = S(A ξ )N + S((A ) ξ )N Hence L = ω * ω * N , [(A ) ξ ∧ N] + e 3 , [(A )ξ ∧ N] dω = H 1 u 1 + H 2 u 2
where

H i = ω * ω * N , [(A ) e i ∧ N] + e 3 , [(A )e i ∧ N] dω, (6a)
3 Using the fact that AS(ξ )ν = S(Aξ )Aν for all ξ , ν ∈ R 3 and A ∈ SO(3).

for i = 1, 2. One may construct the control law

u i (t) = -H i (A(t, •), N(t, •)), i = 1, 2. (6b) obtaining L = -(H 2 1 + H 2 2 ) (7)
The closed loop system ( 5)-(6a)-(6b) with initial condition N 0 (ω) = M 0 (ω) will be called the closed loop auxiliary system.

Main result

Before stating the main result, note that system (5) is parameterized by T in the sense that the adiabatic inputs (4) depends on T , and so is the adiabatic propagator A(t, ω). This implies that the (continuous) solution N(t, ω) of ( 5) depends on the chosen T . The next theorem assures the solution of the problem of Def. 1.

Theorem 2 There exists δ > 0, such that for all ε > 0 and any initial condition M 0 ∈ H 1 ((ω * , ω * ), S 2 ) such that M 0 + e 3 H 1 < δ , then there exist4 T > 0 and ∈ N and a control law

Ω : [0, T ] → R 2 such that M( T, •) + e 3 L ∞ ≤ ε.
Furthermore, for the chosen T and , the control law Ω(t) to be applied to (1) can be obtained in the following way:

Step 1: Determine the adiabatic control

Ω(t) = ( ū(t), v(t)) given in (3).
Step 2: Compute the solution of the closed loop system (5)-(6a)-(6b) with initial condition N 0 = M 0 in the interval [0, T ], obtaining the feedback law (u 1 (t), u 2 (t)) given by (6a)-(6b).

Step 3:

Apply Ω(t) = (u(t), v(t)) = (u 1 (t) + ū(t), u 2 (t) + v(t)) to the system (1) .

COMPUTER SIMULATIONS

One has chosen a priori T = 20, s 0 = 0.1, K = 10, = 16 and T f = T for the new method 5 . As defined in (6a)-(6b), we have chosen unitary gains of the feedback law (we mean, there is no gain multiplying H i of (6b)). For the old method 6 , we have chosen T = 1 and unitary gains as well. We have verified that greater values of T than 1 for the old method are worse, but smaller values of T will not improve the result. Figure 1 shows the simulation results in the Bloch sphere for these choices.

The obtained error of the adiabatic propagator is A(T -, •) -I ≤ 0.0009 obtained in our computer simulation. So ε 2 = A(T -, •) -I ≤ 0.015. It is very small in this case. In the simulations we have found that ε 1 = N( T ) + e 3 L ∞ is more than ten times greater than ε 2 . Hence one will show only the behaviour of the auxiliary state N(t, ω).

In figure 1 one may see the initial condition N 0 = M 0 and the final condition N( T ). One has obtained ε 2 = N(T f ) + e 3 L ∞ = 0.185 with our new method, and ε 2 = 0.58 with our old method. Both methods have considered the same unitary gains multiplying H i (see (6a)-(6b)). The expressions of the feedback of the old method is analogous to (6a)-(6b), with the difference that A(t, ω) is replaced by the matrix exp(Sωe 3 σ (t)), where σ (t) = (-1) E(t/T ) , E(s) is the integer part of s, and σ (0) = 0). Our new method have produced a result that is more than 3 times better than the old method with respect to the final L ∞ norm..

L (t) = 1 2 N(t) + e 3 2 H 1 .
In that figure one shows also the evolution of N(t) + e 3 L ∞ The controls u 1 (t) and u 2 (t) are also depicted in that figure. Figure 3 is a "zoom" of the last one. This allows to see the "microstructure" of the control of the new strategy.

Figure 4 shows the plot of log( N(t) + e 3 2 H 1 ) versus time. The inclination of the curves of log( N(t) + e 3 2 H 1 ) would give a measure of the exponential rate of decaying of N(t) + e 3 2 H 1 = 2L (t). One sees that the convergence is not exponential, the inclination is decreasing as time passes for both methods. However, the inclination is much bigger for the first method in the beginning, and this inclination decreases faster for the old method with respect to the new one. This indicates that the new method seems to be more effective than the old one.

The Figure 5 presents a comparison of the input norms of the old and the new method.

CONCLUSIONS AND FUTURE RE-SEARCH

The main result of this work indicates that the Rabi pulses that are commonly encountered in Nuclear Magnetic Ressonance (NMR) techniques (for instance spinecho pulses) are not a mandatory ingredient for an efficient open loop control law. One might ask if this could imply that one may develop NMR methods with pulses with less intensity than the ones that are found in the present state of the art. This could be an interesting topic of future research, which may lead to produce less "agressive" NMR techniques for medical (and other possible) applications.

ACKNOWLEDGMENTS

The first author is supported by CAPES. The second author is supported by CNPq and FAPESP.

A. Definition of normalized functions a(s), b(s), and k(s)

In this appendix we define the functions s), ( ¯s), and k(s) that are used in the adiabatic control. A computer simulation is presented in order to illustrate the convergence result (Theorem 1). For this, let s 0 ∈ (0, 1/4). Define the function a : [0, 1] → R by (see Figure 6):

a(s) =                -1, if s ∈ [0, s 0 ]; -cos 2π(t-s 0 ) (1-4s 0 ) , if s ∈ (s 0 , 1 2 -s 0 ]; 1, if s ∈ (1/2 -s 0 , 1/2 + s 0 ]; -cos 2π(t-3s 0 ) (1-4s 0 ) , if s ∈ ( 1 2 + s 0 , 1 -s 0 ]; -1, s ∈ (1 -s 0 , 1]. (8a) Define the function b(•) by b(s) = 1 -{ ā(s)} 2 (8b) and k(•) by k(s) = 1, if s ∈ [0, 0.5), -1, if s ∈ [0.5, 1] (8c)
One may extend these functions a, b, k to be 1-periodic functions in a natural way. Figure 6 shows these functions for s 0 = 0.1. A computer simulation of the adiabatic propagator A(t, ω) was done for T = 10, T = 15 and T = 20, with s 0 = 0.1 and K = 10. The values of A(T -, ω) -I as a function of ω is given in Figure 7. The fast convergence of the maximum value of this norm to zero when T → ∞ is easily seen in that figure.

B. Proof Sketches B.1. Auxiliary Results

In order to prove the main result, one proves first the following important results Theorem 3 Fix t 0 = kT for some k ∈ N. Let 7 T 0 > 0 and τ 0 = s 0 T 0 with s 0 ∈ (0, 1/4). There exists δ > 0 with the following property: for all ε > 0, there exist c > 0 such that, for every T ≥ T 0 , and for every initial condition N 0 = N(t 0 ) such that N 0 + e 3 H 1 < δ , and N 0 + e 3 L ∞ ≥ ε, then, one will have L (N(t 0 + τ 0 )) ≤ (L (N 0 )c) for system [START_REF] Beauchard | Semi-gobal weak stabilization of bilinear schrdinger equations[END_REF] in closed loop with the control law (6a)-(6b).

The proof of Theorem 3 is based on the following Lemmas, whose proofs may be found in [START_REF] Beauchard | Stabilization of an arbitrary profile for an ensemble of half-spin systems[END_REF] (particularized for R = I):

Lemma 1 Fix t 0 = kT . There exists δ > 0 such that, if N 0 (ω) + e 3 < δ and the control law (u 1 (t), u 2 (t)) defined by (6a)-(6b) is null on [t 0 ,t 0 + τ 0 ], then N 0 = -e 3 .

The following result is also similar. Since the auxiliary system is T -periodic, we shall state the result for t 0 = 0.

Lemma 2 Consider a sequence of initial conditions

(N n 0 ) n∈N of H 1 ([ω * , ω * ], S 2 ) such that N n 0 N ∞ 0 weakly in H 1 and the associated controls (Ω n ) n∈N , Ω ∞ . Then the solution N n (t, •) N ∞ (t, •) weakly in H 1 and Ω n (t) → Ω ∞ (t) for t ∈ [0, τ 0 ].
Proof. (Of Theorem 3) Let γ = N 0 (ω) + e 3 H 1 < δ , where δ is defined in Lemma 1. Since the auxiliary system is T -periodic, there is no loss of generality in considering t 0 = 0. The proof of this theorem is based on Lemmas 1 and 2. By contradiction, if the result does 7 By construction, for T ≥ T 0 the adiabatic control ( ū, v) is null for t ∈ [0, τ 0 ]). not hold, one may construct a sequence N n 0 , n ∈ N of initial conditions of the auxiliary system with the following properties:

(i) N n 0 + e 3 L ∞ ≥ ε, ∀n ∈ N; (ii) N n 0 + e 3 H 1 ≤ γ, ∀n ∈ N; (iii) τ 0 0 [(u n 1 ) 2 (t) + (u n 2 ) 2 (t)]dt ≤ 1/n, ∀n ∈ N, n > 0. By (ii), passing to a convenient subsequence if neces- sary, one may assume N n 0 N ∞ 0 weakly in H 1 . In par- ticular, N n 0 → N ∞ 0 strongly in the L ∞ norm. Moreover, due to weak convergence, N ∞ 0 +e 3 H 1 ≤ lim n→∞ inf N n 0 + e 3 H 1 ≤ δ [6]
. Now we shall show that the initial condition N ∞ 0 produces null controls for t ∈ [0, τ 0 ], which will be a contradiction with respect to Lemma 1. Now, by the Lemma 2, one has that Ω n (t) = (u n 1 (t), u n 2 (t)) → Ω ∞ (t) where Ω ∞ (t) is the control that is obtained with the initial condition N ∞ 0 . An extra work shows that the controls are of class C 1 , and they are uniformly bounded, as well as their time-derivatives. In particular, the sequence of controls Ω n (t) are uniformly bounded and equicontinuous, and so by Ascoli-Arzela theorem, passing to a subsequence if necessary, Ω n converges to Ω ∞ in C 0 with the sup norm. Assuming that Ω ∞ is not identically null, this gives a contradiction with the fact that the L 2 norm of Ω n tends to zero.

In order to prove the next results, one introduces the Propagator equation of the original system (1): Ḃ(t, ω) = S [u(t)e 1 + v(t)e 2 + ωe 3 ] B(t, ω) [START_REF] Kobzar | Pattern pulses: design of arbitrary excitation profiles as a function of pulse amplitude and offset[END_REF] with initial condition B 0 (ω) = B(0, ω) ∈ H 1 ((ω * , ω * ), SO(3)). By proposition 1 of [START_REF] Beauchard | Stabilization of an arbitrary profile for an ensemble of half-spin systems[END_REF], there exists B 0 ∈ H 1 ((ω * , ω * ), SO(3)) such that B 0 e 3 = M 0 . It is then clear that the solution of the original system (1) is such that M(t, ω) = B(t, ω)e 3 , where B(t, ω) is the solution of (9) and M(t, ω) is the solution of (1) with the same applied input (u(t), v(t)). It is also important to consider the propagator equation that is related to system [START_REF] Beauchard | Semi-gobal weak stabilization of bilinear schrdinger equations[END_REF] given by: Ċ(t, ω) = S A (t, ω)(u 1 (t)e 1 + u 2 (t)e 2 ) C(t, ω) [START_REF] Li | Control of inhomogeneous quantum ensemble[END_REF] One denotes by C(t, ω) the (continuous) solution of [START_REF] Li | Control of inhomogeneous quantum ensemble[END_REF] with the same initial condition C 0 = B 0 and input (u 1 (t), u 2 (t)) such that u(t) = u 1 (t) + ū(t) and v(t) = u 2 (t) + v(t). It is then clear that N(t, ω) = C(t, ω)e 3 , where N(t, ω) is the solution of [START_REF] Beauchard | Semi-gobal weak stabilization of bilinear schrdinger equations[END_REF]. Define the (discontinuous) map B 1 (t, ω) = A(t, ω)C(t, ω), where A(t, ω) is the (T -periodic) adiabatic propagator defined in (3). The following proposition regards the transformations between solutions of (1) and [START_REF] Beauchard | Semi-gobal weak stabilization of bilinear schrdinger equations[END_REF] and between solutions of ( 9) and ( 10):

Figure 2

 2 Figure 2 regards only the new method. It shows the evolution of the Lyapunov functional L (t) = 1 2 N(t) + e 3

Figure 1 .

 1 Figure 1. Results on the Bloch sphere for both the new and the old method.

15

  Lyapunov functional L(t) and N(t) -e3 L ∞

Figure 2 .

 2 Figure 2. Lyapunov functional and inputs for the new method.

15Figure 3 .

 3 Figure 3. Lyapunov functional and inputs for the new method (zoom).

Figure 4 .

 4 Figure 4. Plot of the natural logarithm of the square of the H 1 norm of (N(t) + e 3 ).

3Figure 5 .1 1 aFigure 6 .

 516 Figure 5. Plot of input norms u 1 (t) + u 2 (t) .

Figure 7 .

 7 Figure 7. Plot of the Frobenius norm A(T -, ω) -I as a function of ω with s 0 = 0.1 and K = 10, for T = 10, T = 15 and T = 20.

Recall that, for continuous functions, the L ∞ norm is equivalent to the sup norm.

The proof of Theorem 2 shows that the values of T and of must be big enough.

See Appendix A for the definition of s 0 .

One has simulated the method of[START_REF] Beauchard | Stabilization of an arbitrary profile for an ensemble of half-spin systems[END_REF] considering the trivial target profile identically equal to -e 3 . This is essentially equal to the method of[START_REF] Beauchard | Stabilization for an ensemble of half-spin systems[END_REF].

The first author was supported by CAPES (grant: ) and the second author was supported by FAPESP (grant: ) and CNPq (grant: )

Proposition 1 Assume that M 1 (t, ω) = A(t, ω)N(t, ω). Then, the state N(t, ω) is a solution of [START_REF] Beauchard | Semi-gobal weak stabilization of bilinear schrdinger equations[END_REF] in the interval [kT, (k + 1)T ) with input (u 1 (t), u 2 (t)) if and only if the state M 1 (t, ω) is a solution of [START_REF] Kobzar | Pattern pulses: design of arbitrary excitation profiles as a function of pulse amplitude and offset[END_REF] in the same interval with input u(t) = u 1 (t) + u(t) and v(t) = u 2 (t) + v(t). Similarly, assume that B 1 (t, ω) = A(t, ω)C(t, ω). Then, the propagator C(t, ω) is a solution of [START_REF] Li | Control of inhomogeneous quantum ensemble[END_REF] in the interval [kT, (k + 1)T ) with input (u 1 (t), u 2 (t)) if and only if the propagator B 1 (t, ω) is a solution of [START_REF] Kobzar | Pattern pulses: design of arbitrary excitation profiles as a function of pulse amplitude and offset[END_REF] in the same interval with input u(t) = u 1 (t) + u(t) and v(t) = u 2 (t) + v(t).

In one hand, if one fixes a (continuous) solution B(t, ω) of the propagator equation ( 9), one will obtain a "discontinuous solution" C 1 (t, ω) = A T (t, ω)B(t, ω) of the propagator equation [START_REF] Li | Control of inhomogeneous quantum ensemble[END_REF]. On the other hand, we shall consider (continuous) solutions of the auxiliary system [START_REF] Li | Control of inhomogeneous quantum ensemble[END_REF], that will produce a "discontinuous solution" B 1 (t, ω) = A(t, ω)C(t, ω): Lemma 3 Consider the corresponding (continuous) solution of C(t, ω) of the propagator equation [START_REF] Li | Control of inhomogeneous quantum ensemble[END_REF] that is obtained by the application of an input (u 1 (t), u 2 (t)) in an interval [0, T f ] with initial condition C 0 = B 0 . Let B 1 (•, ω) be the map defined by B 1 (t, ω) = A(t, ω)C(t, ω), where A(t, ω) is the T -periodic map defined in section 2.1. Then the (continuous) solution B(t, ω) of ( 9) that is obtained with with the same initial condition B 0 = C 0 and applying the input (u(t), v(t)) = (u

Proof. See appendix B.2.

B.2. Proof of Lemma 3

The proof of the Lemma 3 is a consequence of the right-invariance of the propagator equation and the fact that the Frobenius norm is invariant by right-and leftmultiplications by orthogonal matrices.

Proposition 2 Fix ω ∈ (ω * , ω * ). The following affirmations are equivalent for the propagator equation ( 9):

• B 1 (t, ω) and B(t, ω) are solutions of (9) in the interval [t 0 ,t 1 ) with the same applied inputs u(t) and v(t).

• There exists a right-translation R(ω)

) and given by R(ω) -I .

The proof of the last proposition is a straightforward consequence of the right-invariance of ( 9), the uniqueness of solutions of ( 9) and the invariance of the Frobenius norm with respect of right-and left-multiplications by matrices in SO( 3). An easy consequence of the last proposition is the following one, which means that, if B(•, ω) is a "continuous solution" of ( 9) and B 1 (•, ω) is a "discontinuous solution" of ( 9), both with the same initial condition B 0 (ω) = B 1 (0, ω) = B(0, ω), then the Frobenius distance between B(•, ω) and B 1 (•, ω) is bounded by the sum of discontinuities of B 1 .

Now note that a (continuous) solution C(t, ω) of ( 10) is transformed into a discontinuous map B 1 (t, ω) = A(ω,t)C(t, ω). From Proposition 1, B 1 (t, ω) is a solution of the propagator equation ( 9) inside the intervals [kt, (k + 1)T ). The "discontinuous" solution B 1 (t, ω) is such that, at each instant t = kT one has a discontinuity of B 1 of a distance B 1 (kT + ) -B 1 (kT -) = A(kT + )C(kT + ) -A(kT -)C(kT -) . By the continuity of the map C(•, ω), then C = C(kT + ) = C(kT -). As the Frobenius distance is invariant by right-and left-multiplications by orthogonal matrices, then A(kT + )C -A(kT -)C = A(kT + ) -A(kT -) = A(T -, ω) -I (recall that A(kT + ) = I). So, after periods of T seconds, Proposition 3 implies that the Frobenius distance B 1 ( T, ω)-B( T, ω) is not greater than the sum of all individual discontinuities. This shows Lemma 3.

B.3. Proof of Theorem 2

Proposition 4 Given ε > 0 and T 0 > 0 there exist some δ > 0 and

Proof. Let c > 0 (that depends on ε) be the constant that is defined in Theorem 3. Let p ∈ N such that L (N 0 )pc ≤ 0. By contradiction, assume that N( T, •) + e 3 L ∞ ≥ ε for all ∈ {0, 1, . . . , p}. Since the Lyapunov functional L (t) is nonincreasing, the repetitive application of Theorem 3 at the instants t = kT for k = 0, 1, . . . , p would give L (N(pT )) ≤ L (N 0 )pc < 0. This is not possible since the Lyapunov functional is always nonnegative. So there must exist some ∈ {0, 1, . . . , p} with the claimed property.

Proposition 5 Fixed k ∈ N, the solution M(t, ω) of ( 1) with initial condition M 0 (ω) is such that

Proof. Consider the original system (1) with initial condition M 0 . Let N(t, ω) be the (continuous) solution of (5) with initial condition N 0 = M 0 . It is clear that

Recall that the "continuous solution" M(t, ω) of the system (1) with initial condition M 0 = N 0 is such that M(t, ω) = B(t, ω)e 3 . Since A(kT, ω) = I, then from Lemma 3 one gets M(kT, ω)

Proof. (of the Theorem 2) Fix ε > 0 and choose ε 1 > 0 and ε 2 > 0 such that ε = ε 1 + ε 2 . It is possible to find ∈ N, T 0 > 0 and a control law Ω T : [0, T f ] → R 2 (depending on T ), with T f = T , in a way that the application of Ω T (t) = (u 1 (t), u 2 (t)) to system (5) furnishes N( , •) + e 3 L ∞ ≤ ε 1 for all T > T 0 . Find T > T 0 big enough (depending on ) such that A(T -, •) -I L ∞ ≤ ε 2 and apply the open loop control (u(t), v(t)) = Ω T (t)+( ū(t), v(t)) to system (1), obtaining M( , •) + e 3 ≤ ε 1 + ε 2 = ε.


Sébastien Cordillet

Nicolas Bideau

Mathieu Ménard

Guillaume Nicolas

Gyro-based lower limb asymmetry during a 4-km time trial on a velodrome among high level female cyclists

Keywords: gyro-based asymmetry, cycling, time-trial, IMU

Introduction

Evaluating pedaling asymmetry in cycling is of major interest since it is often associated to lower efficiency, an increased risk of injury or premature fatigue. To quantify asymmetry, algebraic indices are commonly used to evaluate bilateral differences of discrete variables or through waveform analysis by comparing bilateral patterns [START_REF] Pouliquen | Spatiotemporal analysis of 3D kinematic asymmetry in professional cycling during an incremental test to exhaustion[END_REF]. Whereas some studies quantified bilateral muscle activations of lower limbs during straight and curved sections in track cycling [START_REF] Watanabe | Electromyographic analysis of thigh muscle during track cycling on a velodrome[END_REF], no data were reported concerning in situ kinematic asymmetry adjustments during a cycling time-trial (TT), especially considering women. Thus, the aim of this study was to assess the bilateral asymmetry of bodysegment angular rates during a 4-km TT among female cyclists. In particular, the modifications in asymmetry all along the TT were investigated using lower limb mounted gyroscopes.

Methods

Experimental procedure

After a standardized warm-up, 7 high-level female cyclists, all right-leg-dominant, performed a selfpaced 4-km TT on an outdoor velodrome. Angular rates of bilateral thigh, shank and foot were monitored all along the test using IMUs (Cometa MiniWave Infinity) at a sampling rate of 286 Hz with a full scale of ±2000 deg/sec.

Data analysis

Before the trial, an anatomical calibration procedure was performed in order to align IMU sensor frame to body segment frame following a method previously evaluated [START_REF] Cordillet | Estimation of 3D Knee Joint Angles during Cycling Using Inertial Sensors: Accuracy of a Novel Sensor-to-Segment Calibration Procedure Based on Pedaling Motion[END_REF]. Such a method allows to obtain angular rates in relation to the body segment frames. In order to perform cycle-to-cycle analysis, pedaling cycles were extracted from the medio-lateral angular rate of the thigh mounted IMU by identifying dead centers . Medio-lateral component of angular rate for each body-segment was considered for analysis. In order to compare time delays between right and left patterns of angular rates during pedaling cycles, a normalized cross-correlation technique [START_REF] Gouwanda | Identifying gait asymmetry using gyroscopes-A cross-correlation and Normalized Symmetry Index approach[END_REF] was used to identify similarities between the left and right waveforms as follows:

Where is defined as:

and represent mean values for the right and left instantaneous rotations ( and ). The left angular rates were shifted by 180° to synchronize the initial kinematic cycle at bottom dead center, for both legs.

, expressed as pedal angle lag, indicates the angle lag of the right pattern compared to the left pattern and corresponds to the value for which the correlation coefficient is maximal. A positive indicates that the right pattern is shifted forward in the crank cycle relative to the left pattern. For each body-segment, the maximal crosscorrelation value (

) corresponding to the best overlap between the right and left limb angular rates was considered.

Statistics

A Shapiro-Wilk normality test was employed on and . A non-parametric Friedman test was used to determine the effect of distance (each 1-km interval was compared to the first one). Pairwise Wilcoxon test was employed for the multiple comparisons. Level of significance was set at p<0.05. Table 1 Results of cross-correlations between right and left angular rates waveforms during the 4-km TT.

Results and discussion

Evolutions of angle lags ( ) and maximal correlation coefficient ( ) expressed in degrees. All values were compared with the first kilometer. * denoted significant differences with p value lower than 0.05 (n.s : non significant).

Results of cross-correlation between right and left angular rates waveforms during the 4-km TT are presented in Table 1. Lower correlation was obtained for the thigh. Indeed, while high pattern similarity was observed for the foot and the shank ( ), moderate similarity for thigh was highlighted (

). Furthermore, the similarity between left and right pattern of the shank was not affected by the distance. However, in comparison with the first kilometer, the correlation coefficient for the foot was significantly lowered from the third kilometer. For the thigh, the similarity between bilateral waveforms was significantly lowered during the fourth kilometer. Concerning the angle lag between right and left patterns, a positive sign of was overwhelmingly observed and can be associated to the right leg patterns mostly ahead of time compared to the left leg. More specifically, for the foot, absolute value of increased from the third kilometer. For the shank, significant differences were also observed from the third kilometer, but with a decrease of . Finally, concerning the thigh, a significant difference was quantified from the fourth kilometer. Despite some inter-individual differences, the observed kinematic adjustments, especially the reduction of similarity of side-to-side patterns may reflect a response to fatigue mechanisms that cause movement deterioration during the test [START_REF] Pouliquen | Spatiotemporal analysis of 3D kinematic asymmetry in professional cycling during an incremental test to exhaustion[END_REF].

Conclusions

Considering the angular velocity patterns in the sagittal plane, the asymmetry in the thigh was less affected by the distance and thus by fatigue appearance. The foot demonstrated the highest differences in terms of bilateral pattern similarity but also in terms of angle lag between sides. However previous studies during exhaustive exercises in laboratory conditions demonstrated that the asymmetry of sagittal kinematics was less affected than out of sagittal planes [START_REF] Pouliquen | Spatiotemporal analysis of 3D kinematic asymmetry in professional cycling during an incremental test to exhaustion[END_REF]. Thus, an interesting prospect of this study would be to consider asymmetry analysis in frontal and transverse planes during time trials.

Acknowledgements

cycling team for their participation in the present study.


Maximilien Dupont De Dinechin
email: maximilien.dupont.de.dinechin@ens.fr

Matheus Schuh

Matthieu Moy

Claire Maiza

Scaling Up the Memory Interference Analysis for Hard Real-Time Many-Core Systems

Keywords: response time analysis, algorithm optimization, many-core architectures, real-time systems

proposed an algorithm to compute the impact of interference on memory accesses on the timing of a task graph. It calculates a static, time-triggered schedule, i.e. a release date and a worst-case response time for each task. The task graph is a DAG, typically obtained by compilation of a high-level dataflow language, and the tool assumes a previously determined mapping and execution order. The algorithm is precise, but suffers from a high O n 4 complexity, n being the number of input tasks. Since we target many-core platforms with tens or hundreds of cores, applications likely to exploit the parallelism of these platforms are too large to be handled by this algorithm in reasonable time.

This paper proposes a new algorithm that solves the same problem. Instead of performing global fixed-point iterations on the task graph, we compute the static schedule incrementally, reducing the complexity to O n 2 . Experimental results show a reduction from 535 seconds to 0.90 seconds on a benchmark with 384 tasks, i.e. 593 times faster.

I. INTRODUCTION

Programs running in safety-critical real-time embedded systems must remain predictable in terms of execution time to meet the engineering constraints in their specification. Avionics or autonomous vehicles applications, for example, have analysis and decision making in code heavily coupled to time, so each task in the system must be temporally tightly bounded. Usually, these programs are made of periodic loops that activate tasks and any timing deviation might be propagated causing overlapping issues and even functional failure.

For many reasons (energy, performance, integrity, availability), embedded systems are shifting from single-core to multi/many-cores platforms. A many-core is a type of architecture that typically has hundreds of cores and whose computational power mainly relies on the parallelism level of the programs it runs, in contrast with multi-core processors, where a unique core can be quite powerful on its own. In this work we use the Kalray MPPA-256 [START_REF] Dupont De Dinechin | Time-critical computing on a single-chip massively parallel processor[END_REF] as the evaluation many-core platform, but the algorithm can deal with other arbitration policies.

We are interested in computing a program's global Worst-Case Execution Time (WCET) and analyzing how multiple cores may impact this duration: two tasks running simultaneously in distinct cores cannot be granted access to the memory at the same time, and therefore they slow each other down. Such a slowdown is called interference.

In [START_REF] Graillat | Code Generation for Multi-Core Processor with Hard Real-Time Constraints[END_REF] a framework to develop time-predictable real-time systems for many-core architectures is introduced. It is composed of multiple stages, starting with a dataflow application, which is divided into smaller computational blocks that are compiled into C code, resulting in a DAG of tasks, partially ordered by their dependencies. For each task, the WCET in isolation and number of memory accesses are obtained through a tool such as OTAWA [START_REF] Ballabriga | Otawa: an open toolbox for adaptive wcet analysis[END_REF]. Subsequently the tasks are mapped to cores and ordered. In the final step the release dates and Worst-Case Response Time (WCRT), i.e. WCETs taking interference into account, are computed.

Contribution: This paper presents a new algorithm to compute this last step of the framework in O n 2 time for a program divided into n subtasks. Its implementation is done in Python using the Kalray MPPA-256 as target platform, but conceived with generalization in mind, so new architectures can be integrated. The improvement from previous works [START_REF] Rihani | Many-Core Timing Analysis of Real-Time Systems[END_REF] and [START_REF] Rihani | Response time analysis of synchronous data flow programs on a many-core processor[END_REF] is huge, where an algorithm to solve this problem was showcased, but with a O n 4 time complexity making it intractable for very large task graphs. A long version of this paper with more details on the algorithms and results is available in [START_REF] Dupont De Dinechin | Scaling up the memory interference analysis for hard real-time many-core systems (full version)[END_REF].

Organization: Section II presents the problem, elaborating on the expected input, output and hypotheses assumed. In Section III the original solution from [START_REF] Rihani | Many-Core Timing Analysis of Real-Time Systems[END_REF] is briefly explored before the Section IV where our solution is detailed. To conclude, in Section V a complexity analysis and a performance evaluation of the implementation are made.

II. CONTEXT

A. Interference due to arbitration

Hardware arbiters handle how accesses to a shared resource from different initiators are ordered. Multiple types of arbitration policies exist, serving different purposes, such as timing predictability or throughput. The shift to many-core architectures makes the memory bus arbiter a major influence on the execution time of programs.

A simple, deterministic and starvation-free arbitration policy is the Round-Robin (RR). It gives each initiator an equal grant share in circular order, conditioned to the use of this share. This means that cores access the memory one after another, as long as all of them are requesting to read or write data, otherwise they are skipped.

For instance, assuming a bus size of width 1 word with RR arbitration policy, if three cores have to write 8 words to the memory, the first one writes 1 word, then the second one 1 word, then the third one 1 word and this process is repeated until no core needs to write anymore. In a concrete scenario, the first core to get its access granted suffers no interference, but a very detailed analysis would be needed to know which core is delayed and which one is not. Instead, we consider the worst case in the analysis, i.e. that all cores are delayed. With this policy, all three cores are halted 8+8 times, and assuming that each word access takes 1 cycle, they each receive a total interference of 16 cycles.

B. General description of the problem

To precisely estimate the interference, we need to know the time interval during which memory accesses are performed by each core. For this, we use a time-triggered schedule, where tasks, running on cores, are assigned a release date rel (i.e. the task cannot start before this date even if all its inputs are available earlier), and a WCRT R is computed. As a consequence, we can guarantee the absence of interference between two tasks when their execution interval [rel, rel + R] has no overlap.

Given a Directed Acyclic Graph (DAG) of tasks with dependencies, their mapping and schedule onto cores, their WCETs in isolation, their memory accesses and the bus arbiter description, we need to compute release dates for each of those tasks and the total WCRT of the graph, which accounts for the interference between tasks simultaneously accessing (either reading from or writing to) the shared memory. Additionally, some tasks may have a minimal release date, meaning that they must not be scheduled before that date.

The difficulty in solving this problem is that the release dates and interference values are dependent on each other. This means modifying the release dates of tasks can change how they interfere with each other and a new amount of interference might change the release date of yet to be scheduled tasks. However, once a solution is found, the computed release dates allow to always maintain a precise execution: even if the dependencies of a task are executed faster than their WCETs, the task will not be released before, avoiding unexpected interferences.

Figure 1 shows an example of a task set, its initial schedule (top) and then its final schedule accounting for interference (bottom). The mapping is the following: n 0 → PE0; n 1 , n 2 → PE1; n 3 → PE2 and n 4 → PE3. Their WCETs in isolation are respectively 2, 2, 1, 3 and 2. Moreover, there are minimal release dates defined: t = 0 for n 0 , n 3 ; t = 2 for n 1 and t = 4 for n 2 , n 4 . The amount of memory write accesses can be seen in the DAG on the edges between the nodes. In the timing diagram we can see the interference impact on the release dates and WCRT of the tasks, resulting in a global WCRT of t = 7, instead of t = 6 when the interferences ar ignored.

n 0 n 1 n 2 n 3 n 4 1 1 1 1 1 PE0 PE1 PE2 PE3 n 0 n 1 n 2 n 3 n 4 t = 6 PE0 PE1 PE2 PE3 n 0 I:1 n 1 I:1 n 2 n 3 I:2 n 4 t = 7
In the next section we discuss some non-trivial assumptions that allow us to later develop the algorithm using the basic concepts of the problems described here.

C. Approximations and hypotheses

We assume the following constraints: adding a new task to the program can only increase the interference received by other tasks; and for generality purposes, the interference might be non additive, meaning that the interference between n tasks is not necessarily the sum of the interferences between all pairs. However, some bus arbiters have this additivity property, and exploiting this could simplify and speed up the algorithm for those cases.

Also, we add a conservative hypothesis: when multiple tasks are mapped to the same core, they can be treated as a single big task, summing their WCETs, and memory accesses. This hypothesis empirically outputs less pessimistic release times than a more complex approach consisting in computing all the disjoint sets of tasks interfering with a given one.

III. ORIGINAL ALGORITHM

In [START_REF] Altmeyer | A generic and compositional framework for multicore response time analysis[END_REF], an algorithm is proposed to compute a bound on the delay due to interference for a set of sporadic tasks. It served as an inspiration for the algorithm introduced in [START_REF] Rihani | Response time analysis of synchronous data flow programs on a many-core processor[END_REF], which we improve in this paper. [START_REF] Rihani | Response time analysis of synchronous data flow programs on a many-core processor[END_REF] uses two fixed-point iterations to compute the global response time. The first iteration computes the interference between all tasks with a given set of release dates. The second one adjusts all release dates to respect the dependencies. They are repeated until a stable value for the release dates is found or the deadline is crossed, meaning that the task set is unschedulable.

This algorithm was proved to have a O n 4 complexity [START_REF] Rihani | Many-Core Timing Analysis of Real-Time Systems[END_REF] where n is the number of tasks to schedule, which raises scalability issues. The goal of this work is to reduce this complexity allowing it to be applied to hundreds of tasks. Given the task set and initial release dates, the proposed algorithm works incrementally, by adding tasks one by one to the schedule. The algorithm works with a time cursor t, starting from t = 0 and progressing forward. The tasks are divided into three groups:

• Closed (C): t is after their finish date. These tasks have both their final release date and response times computed. • Alive (A): t is between their release and finish date. These tasks have their final release date, but their response time may be influenced by tasks not yet added to the schedule. • Future: t is before their release date, neither the release date nor their response time is computed.

At each iteration, the cursor t jumps to the nearest end date of the current alive tasks or the minimal release date of future tasks, whichever is smaller. New available tasks, i.e. with all dependencies satisfied, are then scheduled, and the interferences that they add to and receive from the current alive tasks are computed. They cannot interfere with dead tasks, because they do not overlap, and their interferences with future tasks will be computed later in the algorithm, when those are added.

With this approach, when a task is scheduled, its release date is definitively set and, as previously discussed, will not be changed with future tasks. The key idea behind the complexity reduction is that only tasks in A need to be considered in the interference analysis.

Figure 2 captures a snapshot of the algorithm being executed. The vertical dashed red line represents the current time cursor position. Only the solid boxes are considered alive tasks. The dotted boxes on the left are the dead tasks, and the ones on the right are the future tasks.

A. Detailed algorithm

The proposed algorithm is given in Algorithm 1 as pseudo code, and detailed below. The inputs are a task set Γ, their initial release dates Θ and response times R, the number of cores c available in the platform, how the tasks are mapped onto them and a shared memory, that may have distinct arbitrated banks reserved for each core to minimize interference.

In the example from Figure 2, we have Γ = {n 0 , . . . , n 10 }, c = 4 and the mapping is as follows: n 0 , n 1 , n 2 → PE0, n 3 , n 4 → PE1, n 5 , n 6 , n 7 → PE2 and n 8 , n 9 , n 10 → PE3.

The time cursor begins at t = 0, with A, the set of current alive tasks, initially empty. The following steps are then repeated until all the tasks are scheduled (at each step we give the corresponding state in the example from Figure 2): 1) C (closed) is the set of tasks ending at time t. It is simply computed by scanning the current alive tasks, and determining if the end of the task (rel + WCRT) equals t. These tasks are then removed from their reverse dependencies list, allowing tasks depending on these closed ones to start.

Example:

C = n 6 2) A (Alive) ← A -C Example: A = n 0 , n 4 , n 9 3) O (Opening) is the set of tasks opening at time t.
It is computed by scanning the head of the stack of scheduled tasks for each core, and determining whether its dependencies are satisfied and if its minimal release date is smaller than or equal to t.

Example: O = n 7 4) A ← A ∪ O

Example: A = n 0 , n 4 , n 7 , n 9 5) For any destination and source task in A that access the same memory bank, we determine if the source task has already been accounted for in the interferences received by the destination. If not, that interference is recomputed by the bus arbiter function, after adding the source to the list of nodes that the destination interferes with. The interferences are computed separately for each memory bank access from the task τ . The total interference received by the task τ is the sum of those values. 6) t is updated to the minimal value between the next smallest release date of future tasks and the next finish time of alive tasks.

B. Complexity

The size of the set of alive tasks A is bounded by the number of cores. Therefore, we access the linear I BUS function a bounded number of times for each progression of t, and the possible values for t are tasks end dates and their minimal release dates, making it at most 2n. The two nested loops then give an overall complexity, with n tasks, b banks and c cores:

O (c 2 • b • n 2 ).

V. EXPERIMENTAL RESULTS

To compare the old and new algorithm on real world scenarios, we generate random DAGs, using a method proposed by Tobita and Kasahara in [START_REF] Tobita | A standard task graph set for fair evaluation of multiprocessor scheduling algorithms[END_REF], explained and used in the original work by Rihani [6].

This method is called layer-by-layer DAG generation. Tasks on the same layer are assigned to cores in a cyclic way: the n-th task of a layer is assigned to Core(n mod number of cores). Tasks have randomly generated WCET, memory accesses and write operations on tasks of the next layer, respectively between [550, 650], [250, 550] and [0, 100]. Two approaches are used to generate the inputs of the benchmark: fixed N L, in which the number of layers is constant and the layer size increases, and fixed LS, in which the layer size stays the same and it is the number of layers that gets enlarged.

The implementation of the original algorithm is done in C++, while the proposed algorithm is written in Python. This mean that there is an interpreter overhead that negatively impact our results mainly for a small number of tasks.

A linear regression computation on a log×log scale from the benchmark values was done to see if the theoretical complexity goes in hand with the practical outcome. Figure 3 shows the results, where NL4 represents a fixed number of layers of 4, and LS4 a fixed layer size of 4. The bus arbiter function used is the Kalray MPPA-256 RR from [START_REF] Rihani | Many-Core Timing Analysis of Real-Time Systems[END_REF]. The complexity of the proposed algorithm always stay under O n 2 , contrary to Rihani's which exceeds O n 4 and even seem to reach O n 5 in the NL64 and LS64 cases. The benchmark has a timeout that the C++ version easily reaches for more than 256 tasks.

In particular, LS64 and NL64 are the random DAGs configuration values that showcase the biggest difference between the two versions. For LS64 and 256 tasks, the C++ version took 1121.79s and the Python one took mere 4.13s, or 270 times faster. For N64 and 384 tasks, the C++ implementation executed for 535.24s and the Python for only 0.9s, or 593 times faster. VI. CONCLUSION This paper introduces a new algorithm to obtain the release dates and response times of applications in the context of realtime systems implemented on many-core architectures. The revisited version shows a significative complexity improvement to O n 2 , which translates to 593 times faster runtime in our benchmark, in comparison with the original version from [START_REF] Rihani | Response time analysis of synchronous data flow programs on a many-core processor[END_REF]. This allows to accomplish the requirements of modern safetycritical real-time systems, scaling to more than 8000 tasks while maintaining a reasonable execution time.

Figure 1 .

 1 Figure 1. Minimalist program mapped to 4 cores and its timing schedules

Algorithm 1 : 3 C 7 A

 137 Proposed scheduling algorithm Input: Set of release dates Θ = {rel 1 , . . . , reln}, set of response times R = {R 1 , . . . , Rn} of tasks {τ 1 , . . . , τn} Output: schedulable, Θ, R OR unschedulable 1 forall k, S k ← stack of tasks scheduled on core k; A ← ∅; t ← 0; 2 while t < +∞ do ← {τ ∈ A | (τ .rel + τ .WCET + τ .inter) = t}; 4 for τ ∈ C do // τ.rev_deps → tasks that depend on τ 5 for τ in τ.rev_deps do 6 τ .deps.remove (τ ); list of cores c = {0, . . . , c -1} do 10 if S k is not empty then // get top of stack without removing 11 τ _next ← S k .peek(); 12 if τ _next.deps is empty AND min_rel of τ is ≤ t then 13 O ← O∪ {τ }; 14 τ .rel ← t; 15 S k .pop(); // removes top of stack 16 A ← A ∪ O; 17 for τ _dest ∈ A do // task target of mem access 18 for τ _src ∈ A do // task source of access 19 for bank b in banks B do 20 if τ _dest and τ _src both access b then 21 if τ _src not in τ _dest.interfers_with[b] then 22 τ _dest.interfers_with[b].add(τ _src); 23 τ _dest.interferences[b] ← I BUS (τ _dest, τ _dest.interfers_with[b], b); 24 t_next ← +∞; 25 for τ ∈ A do 26 t_next ← min(t_next, τ .rel + τ .WCET + τ .inter); 27 for min_rel in minimal release of future tasks do 28 t_next ← min (t_next, min_rel); 29 t ← t_next;

  For a given processor, b and c are constants, so we may simplify this equation to O n 2 .

Figure 3 .

 3 Figure 3. Benchmark plotted results


Annemiek Van Leendert
email: annemiekvanleendert@visio.org

Michiel Doorman
email: m.doorman@uu.nl

Paul Drijvers
email: p.drijvers@uu.nl

Johan Pel
email: j.pel@erasmusmc.nl

Johannes Van Der Steen
email: j.vandersteen@erasmusmc.nl

Supporting braille readers in reading and comprehending mathematical expressions and equations

Keywords: Braille display, braille reader, algebra, speech synthesizer, TPACK model

Braille readers create an overview of an expression by reading one braille character after another or listening to a voice saying the expression or equation aloud. The specific skills and knowledge that mathematics teachers need to assist this complicated process are often lacking. The purpose of this study was to investigate the effect of an intervention consisting of a course focused on braille display support in combination with Text-To-Speech synthesizer for mathematics teachers on the mathematical performances of braille readers. A quasi-experimental approach was taken to answer this question. Five teachers of an experimental group received the intervention whereas the teacher of the control group did not. Both the experimental group, consisting of 10 braille readers, and the control group, consisting of five braille readers, took a pre-and post-test. The results indicated that there was a very small positive effect of the intervention.

Introduction

Reading and comprehending mathematical expression and equations has always been difficult for students who use braille or braille in combination with speech synthesis as their primary reading medium (hereafter braille readers). Mathematics teachers need to have specific skills and knowledge to assist this complicated process. Braille readers read and comprehend mathematical expressions and equations by touching or hearing in a sequential pattern [START_REF] Millar | Understanding and representing space: Theory and evidence from studies with blind and sighted children[END_REF][START_REF] Millar | Theory, experiment and practical application in research on visual impairment[END_REF]. They need to build an overview by reading one braille character after another or listening to a voice saying the expression or equation aloud. Expressions and equations have to be represented in a linear notation because braille and speech are both linear output modalities [START_REF] Stöger | Accessing and dealing with mathematics as a blind individual: state of the art and challenges[END_REF]. Hence, braille readers have a small perceptual view and have to read expressions in a linear notation. This means that they can't benefit from the layout of an expression that helps the sighted student to understand the structure of an expression at a glance [START_REF] Karshmer | Access to Mathematics by Blind Students[END_REF].

In many countries, braille readers in secondary and higher education use a laptop in mathematics class. Screen Reader software sends the mathematical text on the screen to a refreshable braille display, which is connected to the keyboard of the laptop, or to a Text-To-Speech (TTS) synthesizer. Consider, for example, the next expression:

(1)

In this representation, the expression can't be sent to the braille display or TTS synthesizer. It must first be converted to a linear representation (hereinafter pre-braille notation), e.g. (2x^2 + 1)/3. The notation on the braille display is to a certain degree static. This allows for at least a limited spatial overview [START_REF] Archambault | A software model to support collaborative mathematical work between braille and sighted users[END_REF]. On the other hand, every character on the braille display has to be read individually and sequentially which makes braille reading exhaustive rather than selective [START_REF] Hughes | Movement kinematics of the braille-reading finger[END_REF][START_REF] Hughes | On the nonsmooth, nonconstant velocity of braille reading and reversals[END_REF]. That makes it rather difficult to get an overview over an expression or equation [START_REF] Van Leendert | An exploratory study of reading mathematical expressions by braille readers[END_REF]. Moreover, braille characters have a low redundancy, which means that the characters are hard to distinguish from each other [START_REF] Tobin | Is literacy for blind people under threat? Does braille have a future[END_REF]. This can lead to errors in decoding the braille characters into correct elements of the expression or equation.

In addition to the braille display, the TTS synthesizer is a useful screen reader device for reading mathematical expressions and equations. The Screen Reader that controls the speech and that we used in this study is JAWS, which stands for Job Access With Speech. An advantage of speech over braille is the pace of reading because spoken language can be much quicker comprehended then braille characters. Mathematical expressions and equations can be spelled or read aloud. When an expression is spelled, each element in the expression is approached separately. When an expression is read, a speech dictionary determines the spoken text. This dictionary can be extended with a mathematics vocabulary (also influencing how non-mathematical text is read) and can be further extended with mathematical symbols on so-called verbosity levels. The way in which an expression is spoken aloud depends on the content of the dictionary and the verbosity levels and on the verbosity level at which the braille reader reads. For example, "2(x^2 + 4)" can result in "two x squared plus four" when using a low verbosity level and in "two open bracket x squared plus four closed bracket" when using a high verbosity level. In this case, the meaning is lost when the braille reader uses a low verbosity level. In general, braille readers (should) use a low verbosity level for reading non-mathematical text and a high verbosity level for mathematical text. By default, the settings are not adjusted to mathematical text. In that case, choosing a high verbosity level doesn't help.

In summary, the Screen Reader's software can be adjusted in such a way that (almost) all expressions are spelled and read in mathematical vocabulary. This supports the braille reader in learning the mathematical vocabulary, which is vital for the development of mathematical skills [START_REF] Riccomini | The language of mathematics: The importance of teaching and learning mathematical vocabulary[END_REF]. Moreover, this vocabulary corresponds to the vocabulary used by the mathematics teacher in the classroom. This is particularly useful for braille readers in inclusive classrooms, because individual support is not always available.

Both assistive devices, the braille display as well as the TTS synthesizer, have their strengths and weaknesses. Ideally, using both devices one combines their individual strengths to overcome the weakness of each of the two methods [START_REF] Bernsen | Multimodality Theory[END_REF]. For example, uncertainties in braille can be verified or checked by spelling aloud.

The previous sections demonstrate that reading and comprehending mathematical expressions and equations while using the braille display or the TTS synthesizer is not easy. Teachers often don't know how to guide this complicated process. To map out the qualifications of mathematics teachers, the TPACK model can be used. TPACK includes knowledge of technology (TK), pedagogy (PK) and content (CK), as well as insight into the complex interaction between these knowledge components [START_REF] Mishra | Technological pedagogical content knowledge: A framework for teacher knowledge[END_REF]. The idea behind TPACK is that the pedagogical use of (assistive) technology devices is strongly influenced by the content domains on which these devices are situated [START_REF] Graham | Measuring the TPACK confidence of in-service science teachers[END_REF]. For example, the teacher knowledge required to effectively integrate technology in a mathematics classroom may be very different from that required for a language classroom.

This study aims at improving the integration of the braille display and the TTS synthesizer in the mathematics classroom. We addressed the following research question:

What is the effect of improving the knowledge of braille display and Text-To-Speech synthesizer support of mathematics teachers on braille reader's achievement in mathematics?

We expected that braille readers will better perform in mathematics because the adjusted settings of the Screen Reader software are more appropriate for reading mathematical text and the mathematics teachers have the knowledge to teach the braille readers how to use the braille display and the TTS synthesizer in mathematics lessons.

Methods

Design

In this study, the braille readers' mathematics teachers took part in a professional development (PD) course on integrating the braille display and the TTS synthesizer into mathematics lessons. The first part of the course was aimed at assisting each teacher in adjusting the settings of the Screen Reader software. As a result, the braille readers were able to listen to how expressions and equations were spelled and read aloud in mathematical vocabulary. The second part of the course was aimed at improving the knowledge to teach the braille readers how to use the braille display and the TTS synthesizer in mathematics lessons. In order to answer the proposed research question, a quasiexperimental approach was taken, using a pre-and post-test. The experimental group consisted of 10 braille readers whose teachers took part in the PD course. The control group consisted of five braille readers whose teacher did not participate in this course.

Context and participants

In the Netherlands, more than 50% of all braille readers in secondary education go to special schools for students with a visual impairment. Organization A and B provide education for these students. The first author was affiliated with organization A. Therefore, five mathematics teachers from two different schools of organization A participated in the PD course, while the teacher who worked at the school of organization B did not. None of the teachers ever received any training in teaching mathematics to braille readers. The mathematics teachers who worked at organization A were qualified teachers, but they were not, with the exception of one, qualified to teach mathematics. In contrast, the mathematics teacher who worked at organization B was a qualified mathematics teacher. The braille readers of organization A were assigned to the experimental group and the braille readers of organization B were assigned to the control group. This resulted in an experimental group of 10 and a control group of five braille readers. The three schools were very similar in terms of demographic characteristics of students and the degree of educational performance. The braille readers were all in seventh to twelfth grade in secondary education. They differed in their mathematical skills and in their ability to use the braille display and the TTS synthesizer.

Intervention: a professional development course

Part of the intervention consisted of the adjustments of the settings of the Screen Reader software. The control group worked with the "old" settings, and (often) used one verbosity level when reading. As a result, too few symbols and punctuation marks were read aloud in mathematical text, or too many punctuation marks were read aloud in non-mathematical text. Moreover, (almost) no mathematical vocabulary was used. For the experimental group, the verbosity settings of the Screen Reader were adjusted. The speech dictionary was extended with a mathematical vocabulary and all punctuation marks and symbols, used in pre-braille notation, could be spelled aloud in mathematical vocabulary. The braille readers could choose a verbosity level with which they could read the whole expression aloud in mathematical vocabulary, without missing elements.

The face-to-face PD course consisted of four sessions of three hours each, which lasted over a period of four months. The design of the course was based on TPACK. Much attention has been paid on how to integrate the braille display and the TTS synthesizer into the mathematics lessons. For the first homework task, the mathematics teachers were asked to support a braille reader in obtaining an overview over an expression or equation. For the second homework task, they were asked to write down their experiences with the adjustments of the Screen Reader software and to discuss this with their braille readers. For the third task, they were asked to do an activity to provoke collaborative work and mathematical communication between students.

Pre-and post-test

We decided that we could use identical pre-and post-tests because the time interval was more than four months. The braille readers were asked to answer the two tasks orally, because that would save time and be less tiring for them. The tasks were not too complex but required careful reading skills due to the use of various operations and brackets. In the first task, "synthetic speech comprehension", the braille readers had to select information from expressions or equations that were spoken aloud. The items could be spoken aloud several times at the request of the braille reader. The items were: a) 1/7 + 7 -2/5 + 8 = .. What are the fractions in this expression? b) 4 -(5 + -(6 + 3)) = .. Where do you start calculating? c) 4 * (.. -5) = 8 Solve this equation.

In the second task, "mathematical braille reading skills", the braille readers were asked to read the mathematical text on the braille display and verbalize this text in mathematical vocabulary. With item b we expected to be able to investigate whether they were supported by context. The items were:

a) y = 2 1/2 * 3
b) The volume is 12 m^3 c) y = sqrt(2/(x + 3)^2)

(2 1/2 is a mixed number, a combination of a whole number and a proper fraction). Before the start of the pre-and post-tests, the braille readers were interviewed about their visual impairment, their assistive devices and the support they receive in mathematics lessons. This information helped to interpret the results of the tests.

Procedure

The research followed this order: pre-interview, pretest, intervention, post-interview, posttest, all within a maximum of five months. The interviews and the pre-and post-tests took place at the schools of the braille readers. All sessions, consisting of an interview and a pre-or post-test, were scheduled for 25 minutes but sometimes lasted longer due to technological problems with the braille readers' laptops. Within two weeks of completing the pretests, the Screen Reader software of the experimental group was adjusted and the mathematics teachers' development course started. Within a month after the end of the course, the braille readers did the posttest. The control group received the same interviews and pre-and post-tests in the same period but without intervention.

Data collection and analysis

For investigating the research question data was collected during the pre-and post-test through audio and video recordings. For "synthetic speech comprehension", data were collected on whether the braille readers gave a correct answer and on how much time they needed to give this answer. For each group, the average time was calculated. For "mathematical braille reading skills", data were collected on whether the braille readers could verbalize the expression, how much time they needed to give a correct answer and on the kind of errors they made. For each group, the average time was calculated.

Results

The results in Table 1 show that for synthetic speech comprehension, the pretest percentage of correct answers was 57% for the experimental and 47% for the control group. For the posttest, the percentage of correct answers was (again) 57% for the experimental and 80% for the control group.

The experimental group needed in the pretest, on average 46.1 seconds (SD = 27.0) and in the posttest, on average, 34.1 seconds (SD = 16.0) to give the correct answers. The control group needed in the pretest, on average, 45.4 seconds (SD = 25.1) and in the posttest, on average, 32.4 seconds (SD = 14.3) to give the correct answers. For both groups, the standard deviations were high.

For mathematical braille reading skills, the pretest percentage of correct answers was 20% for the experimental and 53% for the control group. In the posttest, the percentage of correct answers was 37% for the experimental and 60% for the control group. The experimental group needed in the pretest, on average, 18.2 seconds (SD = 26.8) and in the posttest, on average, 15.3 seconds (SD = 8.7) to give the correct answers. The control group needed in the pretest, on average, 9.0 seconds (SD = 3.1) and, in the posttest, on average 9.8 seconds (SD = 8.4) to give the correct answers. Also for this task, the standard deviations were high. The errors were categorized. E(d) is an error made due to difficulties with decoding the braille characters. This could be ascribed to difficulties with recognizing the location of the raised dots or to difficulties with decoding the braille characters into correct elements of a mathematical expression (e.g., "x" "^" "2"). E(v) is an error due to difficulties in recognizing and verbalizing the whole expression (e.g. "x squared") and E(b) is an error due to malfunctioning of one of the technology devices. For both groups and for all tests, almost all errors were decoding errors. For both groups, most decoding errors occurred at decoding symmetric or translated characters (e.g. and ) and mixing up six dot braille with eight dot braille. The E(b) error occurred only, two times, in the pretest of the experimental group. Finally, we illustrate the decoding practices of two braille readers of the experimental group reading "The volume is 12 m^3" on the braille display (see Table 2). In this example, I. tried to read 12 m^3 as one number "2393". He did not recognize the "^" character in braille and stopped. F. also started with reading 12 m as one number, but corrected the error. Both braille readers mixed up "m" with "3". and spent a lot of time reading "12 m^3" compared to reading "The volume is".

Conclusions and discussion

This study aims at improving the integration of the braille display and the TTS synthesizer in the mathematics classroom. For "synthetic speech comprehension", the percentage of correct answers for the experimental group (57%) did not change during the intervention. This means that, for this the task, there was no evidence that the intervention was successful. For the control group, this percentage was 47% in the pretest and 80% in the posttest. It is remarkable that this group achieved much better on the posttest. For "mathematical braille reading skills", the experimental group improved more than the control group. It is remarkable that the pretest percentage of correct answers of the experimental group (20%) was much lower than that of the control group (53%). In the posttest, the percentage was 37% for the experimental and 60% for the control group. For both groups, most errors were errors in decoding the braille characters showing that braille readers already stumble in the first phase of the solving process. Table 2 shows that braille reader I. and F. needed more time and made more errors when reading mathematical text ("12 m^3") compared to reading non-mathematical text ("The volume is"). Both braille readers mixed the symmetric characters "3" and "m", which is a common mistake [START_REF] Tobin | Is literacy for blind people under threat? Does braille have a future[END_REF]. In summary, the results of the tasks on mathematical braille reading show a very small positive effect of the intervention. We expected a greater effect. Finally, the high values for standard deviation, in all tests, showed that the individual differences were large.

This study had some limitations. A first limitation was the small number of braille readers that participated in this study, especially the number of participants in the control group. A second limitation was the limited information we had of the braille readers about the frequency that they used the adjustments of the Screen Reader software in the period between the two tests. From a teacher's point of view, we don't know, exactly, if and how they changed their daily practice. This feedback was asked during the second interview, yet the braille readers found it difficult to reflect on this.

Moreover, there were differences between the experimental and control group that may have affected the results. Firstly, the braille readers of the control group were obliged to use the braille display, whether or not in combination with the TTS synthesizer, during their mathematics lessons. This was not obligated for the braille readers of the experimental group. Secondly, the mathematics teacher of the braille readers of the control group was a qualified mathematics teacher. This was not the case for the mathematics teachers, with the exception of one, of the experimental group.

We did not anticipate that braille readers would make so many decoding errors while reading mathematical text. Overall, the control group performed better on both tasks than the experimental group. This may be due to the strict rules for the use of the braille display and the differences in background of the mathematics teachers. This could be examined in more detail in future studies.

Overall, this study adds to the small number of studies into ways to support braille readers in mathematics (e.g. [START_REF] Figueiras | A touch of mathematics: coming to our senses by observing the visually impaired[END_REF]. Findings from this type of research should enable teachers to better support braille readers in doing mathematics.

Table 1 : Results from pretest and posttest for synthetic speech comprehension and mathematical braille reading skills

 1

			Pretest	Posttest
	Condition	Number	Percentage of	Average time on	Percentage of	Average time on
		of	correct	correct answers	correct	correct answers
		students	answers	(s.)	answers	(s.)
	Synthetic	E (10)	57%	46.1 (SD = 27.0)	57%	34.1 (SD = 16.0)
	speech comprehension	C (5)	47%	45.4 (SD = 25.1)	80%	32.4 (SD = 14.3)
	Mathematical	E (10)	20%	18.2 (SD = 26.8)	37%	15.3 (SD = 8.7)
	braille reading skills	C (5)	53%	9.0 (SD = 3.1 )	60%	9.8 (SD = 8.4)
	E denotes experimental and C control group

Table 2 : Results of two individual braille readers on mathematical braille reading skills

 2

		Braille reader I. (pretest)	Braille reader F. (posttest)
	Time	Utterances	Time	Utterances
	(s)		(s)
	00-04	The volume is	00 -01	The volume is
	04 -17 17 -22	Two thousand three hundred ninety-three Wait two thousand What should this be? The volume is	01 -19	Eh hundred three and twenty No no twelve I think this is cubic 12 cubic meters
	22 -40	one, two three no wait a two m s three I
		do not know if that is an s. It looks a bit
		strange ... {stop}

Acknowledgments

This work is part of the research program 'Improving the mathematical abilities of Braille-dependent students' (023004048), which is financed by the Netherlands Organization for Scientific Research (NWO). This study was approved by the medical ethical committee of the Erasmus Medical Centre (MEC-2012-097 and MEC-2012-524) and adheres to the tenets of the Declaration of Helsinki (2013) for research involving human subjects. Informed consent was obtained from the participants.


Denis Laxalde

Jean-Pierre Lombard

Fabrice Thouverez

Dynamics of Multistage Bladed Disks Systems

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

The mechanical design of gas turbine rotors is traditionally fulfilled through stage by stage analysis. However, in modern engines, dynamical phenomena involving multiple stages (adjoining or not) of bladed disks tend to appear and can then be of primary interest for gas turbine professionals. On the other hand, the traditional (single stage) analysis of bladed disks dynamics is usually performed using the cyclic symmetry reduction which allows to analyze a complete rotationally periodic structure by considering only one of its elementary sector. Examples in the literature are numerous on this subject, from the early developments of cyclic symmetry method [START_REF] Thomas | Dynamics of Rotationally Periodic Structures[END_REF]; [START_REF] Wildheim | Excitation of Rotationally Circumferentially Periodic Structures[END_REF] to specific researches on bladed disks dynamics [START_REF] Srinivasan | Vibrations of bladed disks assemblies -A selected survey[END_REF]; [START_REF] Slater | Forced response of bladed disks assemblies -A survey[END_REF]. Withal, in such analysis, the interstage coupling is necessarily neglected, its eventual effects can not be accounted for and this can lead to significant deviations from the real (multi-stage) system. For example, Bladh et al. [START_REF] Bladh | Effects of Multistage Coupling and Disk Flexibility on Mistuned Bladed Disk Dynamics[END_REF] have show that the interstage coupling can significantly affect the dynamics of the multistage assembly and that the traditional single stage analysis can, in some cases, leads to an underestimation of vibratory levels.

Few strategies exist to answer the challenge of multistage modeling and a major issue is that, since a multistage assembly of cyclic structures is a priori no longer cyclic, the cyclic symmetry reduction cannot be used straightforward. Some [START_REF] Rzadkowski | Forced Vibration of Several Bladed Discs on the Shaft[END_REF] have used full 360 degrees models to study the (free and forced) dynamics of multi-stage systems. However, as such techniques often yield prohibitive computational costs. Component mode synthesis methods (with multilevel reductions) has also been used [START_REF] Bladh | Reduced Order Modeling Techniques for Dynamic Analysis of Mistuned Multi-Stage Turbomachinery Rotors[END_REF]; [START_REF] Song | Multi-stage modeling of turbine engine rotor vibration[END_REF] but these strategies can be somewhat difficult to use in a design process.

In the present paper, we put forward a new strategy to study the dynamics of multistage assemblies using a cyclic symmetry reduction. This method handles the sector finite element modeling directly and involves a realistic interstage coupling based on a cyclic recombination of the interface on the chosen nodal diameter. Following a brief review of the cyclic symmetry method (section 2), the proposed method will be detailed in section 3. Then, some numerical examples and results will be presented in section 4 to validate the method and highlight the importance of such multi-stage calculations in the design of turbomachinery components.

Cyclic symmetry for single stages

In this first section, a brief review of the classical cyclic symmetry modeling is done which mainly aims at defining some notations. A cyclic structure composed of N sectors is considered. Its displacement vector can be split as,

x = x 0 , x 1 , . . . , x N -1 T (1)
where x j is the displacement vector of the j-th sector.

As a consequence of the rotational periodicity, the mode shapes are such that each sector has the same deflection and a constant phase difference from the adjacent sectors. This inter-sector phase angle is β n = 2πn/N and the harmonic index n = 0, 1, . . ., N -1 is called the nodal diameter number.

As a result, any dynamical problem on the complete cyclic structure,

M ẍ + C ẋ + Kx = f e (2)
can be transform (exactly) into reduced problems on the reference sector (corresponding to each nodal diameter number n)

M n ẍn + C n ẋn + K n xn = f e,n (3)
The displacement vector, in the cyclic base, xn is obtained from the reference sector's displacement vector in which the left inter-sector degrees-of-freedom, l x 0 , have been eliminated using the relation to the right inter-sector degrees-of-freedom, r x 0 , for each n nodal diameter,

l x 0 = e -2iπn N r
x 0 (4)

The reduced matrices M n , C n and K n are obtained by applied this periodicity conditions on the reference sector's structural matrices.

As the considered structures have symmetric structural matrices, xn and xN-n contains the same relevant information and thus among the N possible values of the nodal diameter index n only N s of these are to be calculated,

N s = N/2 + 1 if N is even (N + 1)/2 if N is odd (5)
Each reduced problem (3) is solved in cyclic coordinates xn (for the reference sector) and the transformation in physical coordinates (complete structure) is given by: where ⊗ is the Kronecker product.

For further readings, see references [START_REF] Thomas | Dynamics of Rotationally Periodic Structures[END_REF]; [START_REF] Wildheim | Excitation of Rotationally Circumferentially Periodic Structures[END_REF]; [START_REF] Henry | Calcul des fréquences et modes des structures circulaires[END_REF].

A multistage cyclic symmetry method

We suggest a new strategy to analyze the dynamics of multistage cyclic structures. This method was originally proposed by the authors in Laxalde et al. (2007a). It is called multistage cyclic symmetry since the cyclic symmetry strategy is used to model each stage. The interstage coupling is ensured through a cyclic expansion of the interface degrees of freedom of each stages. The method will be detailed based on the example of two stages respectively composed of N 1 and N 2 sectors (N 1 < N 2 ). This method is applicable for any dynamical problems in free and forced responses. In this paper, we focused on linear applications for modal analysis and forced steady-state response. Extension to non-linear steady-state response calculations can easily be considered by combining cyclic symmetry approach with Fourier methods [START_REF] Petrov | A method for use of cyclic symmetry properties in analysis of nonlinear multiharmonic vibrations of bladed disks[END_REF]; [START_REF] Laxalde | Vibration control for integrally bladed disks using friction ring dampers[END_REF].

Equations of the motion. We consider the general dynamical problem involving several coupled stages, each having a second-order equation of the motion of the form,

M (s) ẍ(s) + C (s) ẋ(s) + K (s) x (s) = f (s) e (s = 1 or 2) (7)
in which, M (s) , C (s) and K (s) are mass, damping and stiffness matrices, x (s) is the displacement vector and f

(s)
e is the external forces on the full stage s. Since each stage is rotationally periodic, it can be represented by its elementary sector in its own cyclic base by the following uncoupled nodal diameters equations,

M (s) n ẍ(s) n + C (s) n ẋ(s) n + K (s) n x(s) n = f (s) e,n , for n ∈ [0 . . . N (s) s -1] (8)
The proposed method aims at coupling the models of each stage in their own cyclic components.

Interstage coupling. At the interface between stages 1 and 2, the coupling is achieved by enforcing continuity of displacements (or other variables) on the interstage boundary in the common physical base,

∆ b x = A b x (1) -b x (2) = 0 (9)
where the left superscript b . refers the restrictions to the interstage boundary degrees-of-freedom and A is a constraint matrix which makes the two interstage meshes compatible (multi-point constraints). Equation. ( 9) can be rewritten, using each stage's cyclic base, as,

A N 1 -1 n=0 e N 1 ,n ⊗ b x(1) i,n + N 2 -1 p=0 e N 2 ,p ⊗ b x(2) i,p = 0 (10)
Then, we can project equation ( 10) on any nodal diameter vector of stage

2 cyclic base (p ∈ [0 . . . N 2 -1]) and obtain, e * N 2 ,p ⊗ I b 2 A N 1 -1 n=0 e N 1 ,n ⊗ b x(1) i,n + b x(2) i,p = 0 (11)
where the identity mtrices I bs (s = 1, 2) ensure the compatibility of the interstage interface.

In equation ( 11), the coupling terms between nodal diameters n and p of each stage's base are

C(p, n) = e * N 2 ,p ⊗ I b 2 A e N 1 ,n ⊗ I b 1 (12)
The cornerstone of the proposed method is to neglect the coupling between "distinct" nodal diameters of each stage's base. We can, here, define the correspondence of nodal diameters indexes between the two stage's bases (which may not be obvious since the stages have different number of sectors). For any nodal diameter index p ∈ [0 . . . N 2 ] of stage 2, we can define the corresponding "equivalent" nodal diameter index of stage 1 n(p) as,

n(p) =      p if p ≤ N 1 N 1 - p (mod N 1 ) if N 1 is odd p (mod N 1 ) if N 1 is even if p > N 1 (13)
This correspond to the aliasing of the discrete Fourier transform. The "nodal diameter equivalence" is illustrated by figures (1), with the even (fig. 1

* N 2 ,p ⊗ I b 2 A e N 1 ,n(p) ⊗ I b 1 b x(1) i,n(p) + b x(2) i,p = 0 (14)
This assumption will next be referred as the one-to-one nodal diameter coupling assumption.

Multistage cyclic symmetry reduction. The one-to-one coupling assumption makes possible the multi-stage coupling in the cyclic bases for each nodal diameters separately. We can append two stages models in their own cyclic base and according to the previously defined nodal diameter correspondence relation ( 13) as,

M (1) n(p) 0 0 M (2) p ẍ(1) n(p) ẍ(2) p + C (1) n(p) 0 0 C (2) p ẋ(1) n(p) ẋ(2) p + K (1) n(p) 0 0 K (2) p
x( 1)

n(p)
x( 2)

p = f (1) e,n(p) f (2) e,p (15)
Next, if each stage displacements vector is split as,

x(1) n(p) = b x(1) n(p) , i x(1) n(p) T and x(2) p = b x(2) p , i x(2) p T
Eq. ( 14) can be used to eliminate the cyclic components of the interstage boundary of stage 2,

x( 1)

n(p)
x( 2)

p = T p x(1,2) p (16)
with,

T p =     I b 1 0 0 0 I i 1 0 B p 0 0 0 0 I i 2     , B p = e * N 2 ,p ⊗ I b 2 A e N 1 ,n(p) ⊗ I b 1
and,

x(1,2)

p = b x(1) n(p) i x(2) n(p) i x(2) p T
The final multi-stage one-to-one nodal diameter dynamical problem can be defined as:

M (1,2) p ẍ(1,2) p + C (1,2) p ẋ(1,2) p + K (1,2) p x(1,2) p = f (1,2) e,p (17)
with,

M (1,2) p = T T p M (1) n(p) 0 0 M (2) p T p C (1,2) p = T T p C (1) n(p) 0 0 C (2) p T p K (1,2) p = T T p K (1) n(p) 0 0 K (2) p T p
This method allows the study of nodal diameter modes separately as for individual stages, each being modeled with its own elementary sector. The size of problems ( 15) or ( 17) is close to the sum of the two individual stages problem sizes. The accuracy of this method in modal analysis and particularly the validity of the one-to-one nodal diameter coupling will be demonstrated by numerical examples in section 4.1.

Superposition. The multistage cyclic symmetry formalism allows to isolate the nodal diameters contributions of the response. However, in some cases and particularly in forced response, the dynamics of the multi-stage system can feature simultaneous contributions of different nodal diameters. To handle with such problems, two strategies can be used. The first is to use the "superposition principle" as in typical single stage cyclic symmetry analysis; that is to say, the external forcing distribution on each stage is first decomposed into elementary cyclic forces corresponding to nodal diameter contributions of each stage and then each uncoupled, nodal diameter projected, equation can be solved individually. The global response is then obtained by summing all contributions. The second is to consider all harmonic components of the excitation (or response) in the projection step of the multi-stage cyclic symmetry method. Example given, if two harmonics (or nodal diameter numbers) p and q are considered, the projection equation ( 14) becomes, B p,q

x(1) i,n(p)

x(1)

i,n(q) + x(2) i,p x(2) i,q = 0 (19a) with B p,q = e N 2 ,p , e N 2 ,q * ⊗ I b 2 A e N 1 ,n(p) , e N 1 ,n(q) ⊗ I b 1 (19b)
4 Numerical results

The next of the paper is dedicated to numerical applications. First, a validation example will be used to evaluate the performances of the proposed method in modal analysis. Then, an industrial application based on two HP compressor rotating stages will be presented in both modal and forced responses. The importance of multi-stage dynamic analysis will also be discussed.

Validation example

In this part we evaluate the performance of the Multistage cyclic symmetry method on a numerical example. We consider two stages of cyclic structures with simple geometry; the first being composed of 24 elementary sectors and the second of 45 elementary sectors. The figure 2 shows the finite elements meshes of the two elementary sectors and of the whole structure. The results on the Multistage cyclic symmetry method are compared to a full 360 degrees analysis taken as a reference.

Correlation on modal data Figure 3 gathers the eigenfrequencies calculated by the two methods plotted versus their nodal diameter index. The correlation appears to be correct. This attests for the quality of the prediction of the multistage cyclic symmetry method and the validity of the one-to-one nodal diameter projection assumption. Next, we show and compare (in a qualitative way) the modal deformed shapes. Figures 4 and5 show results from (a) the reference analysis and (b) the Multistage cyclic symmetry for various nodal diameter numbers. In these examples we can distinguish some modes localized in only one stage (figure 5) and some others with a clear multistage behavior (figure 4). Here again, the correlation is correct.

Correlation on frequency responses

We now investigated the forced response of the two stages model of figure 2. We consider two distributions of forces on each blade of each stage with a spatial periodicity of n (n nodal diameters), such that the force of the j-th blade of stage s is, For a wave number of n = 2, the frequency responses are displayed in figures 6 for the first blade of each stage. The correlation appears to be globally correct. The resonances (frequency and levels) are correctly represented by the multi-stage cyclic symmetry. There only appears small shifts on the antiresonances between the two calculations. This might be related to the one-to-one approximation of the multi-stage cyclic symmetry method described in section 3. This is quite similar to a truncation effect in modal analysis. However, the global dynamics remains well approximated.

f (s) j = F 0 cos 2nπ N s (j -1) cos ωt, s = 1 or 2 (20)

Industrial case study

Next, an industrial application of the multi-stage cyclic symmetry method is presented. Figures 7 show the finite element models (sectors) of two HP compressor blisks which respectively have 36 and 60 sectors. The computation times of a traditional multi-stage analysis (with 360 degrees models) are clearly prohibitive and thus useless during a realistic design process; in contrast, multi-stage cyclic symmetry analysis can provide satisfying results with reasonable computational times.

Modal analysis. For this representative model, we found interesting to show the influence of the interstage coupling by comparing the results of single stage analysis to the multistage case on a modal analysis. For the single stage analysis we applied a fixed boundary condition on the interstage interface. Figure 8 gathers the (normalized) eigenfrequencies of these three calculations. The influence of the interstage coupling is noticeable and new modes appear that could not be predicted by a single stage analysis. This is particularly observable in small nodal diameter numbers. In effect, in these cases, the blade/disk coupling is generally more important than for higher nodal diameter numbers where the disk's participation is usually smaller. In addition, we provide some example of modal deformed shapes from the multistage computation in figure 9. Some modes show a weak interstage coupling, such as the ones represented in figures 9(a) and 9(b), and feature nearly a single stage behavior. On the other hand, figure 9(c) displays a strongly (interstage) coupled behavior and can obviously not be represented in traditional single stage analysis. This comparison clearly shows the influence of the multistage coupling in the design of multi-stage blade disks systems. (with fixed inter-stage interface) to the two stages coupled. These response were calculated using a modal superposition and a nodal diameter superposition as explained in section 3.

The first remark would be that the multi-stage responses have a higher frequency content than the single-stage responses and that some peaks appears due to the interstage coupling. This is a major remark for the designer. In effect, as the bladed disks are usually designed by avoiding any coincidence of the modal response to the engine excitation, the possible transfer of vibratory energy between stages (adjacent or not) can, with the present method, be predicted. An other important point relates to the vibration levels which significantly differ between the two computations; in effect, even when the peaks location (frequency) seems to be correctly evaluated by the single stage analysis, the resonant level are generally lower in the multi-stage case. Because of the inter-stage coupling, the vibratory energy seems to spread over the two stages and this tends to reduce the vibratory levels in individual stages.

Conclusions

A new strategy to study the dynamics of multistage cyclic structures, such as turbine engine bladed disks, has been proposed. It is based on a cyclic symmetry modeling of individual stages (in their respective cyclic base) and an accurate interstage coupling in the common physical base. This method appears quite simple to carry out since the sector models of each stage are handled directly and the coupling involves few additional difficulties. The results from numerical simulations show a good accuracy of this method. Also, an industrial application highlight the importance of such calculations for the designer.

Acknowledgments

Thanks go to Snecma for its technical and financial support. This work takes place in the framework of the MAIA mechanical research and technology program sponsored by CNRS, ONERA and SAFRAN Group.

x

  = e n ⊗x n with e n = 1 e

  (a)) and odd (fig. 1(b)) cases. Finally, the N1 = 13 and N2 = 20

Figure 1 .

 1 Figure 1. "Nodal diameter equivalence" rule

Figure 2 .

 2 Figure 2. Finite element model of two cyclic structures -N 1 = 24, N 2 = 45

Figure 3 .

 3 Figure 3. Eigenfrequency correlation, reference / Multistage cyclic symmetry + .

Figure 4 .Figure 5 .

 45 Figure 4. One nodal diameter mode at 4824Hz

Figure 6 .

 6 Figure 6. FRF correlation; (---) reference, (--) multi-stage cyclic symmetry

Figure 8 .

 8 Figure 8. Influence of interstage coupling, Eigenfrequency/nodal diameter map (from 0 to 10 ND), stage 1 alone ×, stage 2 alone +, coupled stages .

Figure 9 .Figure 10 .

 910 Figure 9. Examples of modal deformed shapes