Magalie Thomassin
email: magalie.thomassin@laps.ims-bordeaux.fr

Rachid Malti
email: rachid.malti@laps.ims-bordeaux.fr

MULTIVARIABLE IDENTIFICATION OF CONTINUOUS-TIME FRACTIONAL SYSTEM

This paper presents two subspace-based methods, from the MOESP (MIMO output-error state space) family, for state-space identification of continuous-time fractional commensurate models from sampled input-output data. The methodology used in this paper involves a continuous-time fractional operator allowing to reformulate the problem so that the state-space matrices can be estimated with conventional discrete-time subspace techniques based on QR and singular value decompositions. The first method is a deterministic one whereas the second approach takes place in a stochastic context. The performance of both methods is demonstrated using Monte Carlo simulations at various signalto-noise ratios. The deterministic method leads, as expected, to biased estimates. This bias is removed in the stochastic method by the use of an instrumental variable. As compared to rational systems, the commensurate differentiation order must be estimated besides the state-space matrices which is done using nonlinear programming. This is the first work developed for multiinput multi-output system identification using fractional models.

INTRODUCTION

.

. They are based on the minimization of an output error criterion by nonlinear programming techniques. These methods are well suited for single-input singleoutput (SISO) systems, and are generally difficult to apply in the multi-input multi-output (MIMO) case because the number of parameters to estimate becomes large. Here, subspace methods are proposed to estimate the matrices of the continuous-time 1

fractional state-space representation. It is an extension of the methods presented in the literature for rational (thus non fractional) systems [START_REF] Haverkamp | Identification of continuous-time MIMO state space models from sampled data, in the presence of process and measurement noise[END_REF][START_REF] Johansson | Stochastic theory of continuous-time state-space identification[END_REF][START_REF] Johansson | Stochastic theory of continuous-time state-space identification[END_REF] to the fractional case. Other subspace techniques for identifying continuous-time rational models can be found in [START_REF] Bastogne | A PMFbased subspace method for continuous-time model identification. Application to a multivariable winding process[END_REF][START_REF] Ohsumi | Subspace identification for continuous-time stochastic systems via distribution-based approach[END_REF][START_REF] Li | Subspace identification of continuous time models for process fault detection and isolation[END_REF][START_REF] Mercère | Subspace-based methods for continuous-time model identification of MIMO systems from filtered sampled data[END_REF]. So, the proposed method inherits the advantages of subspace methods which stem from the reliability of numerical algorithms using the QR and the singular value decompositions [START_REF] Katayama | Subspace methods for system identification[END_REF]. Thus, it does not involve nonlinear optimization to obtain state-space matrices. In addition, no canonical form (such as modal or companion realizations) of the statespace representation is required. Finally, the proposed subspace algorithms can be applied to the identification of both SISO and MIMO fractional systems. As will be seen later, the state-space representation of a fractional commensurate system involves an additional parameter which is the commensurate order. It is the only parameter computed by minimizing an output error criterion with a nonlinear optimization technique.

In section 2, some recalls about fractional systems are presented. Section 3 presents the methods proposed to estimate the matrices of the continuous-time fractional state-space representation, followed by section 4 devoted to the estimation of the fractional commensurate order using a nonlinear optimization technique. Finally, simulation examples are given in section 5. Monte Carlo simulations are made to show the estimator statistical properties.

FRACTIONAL SYSTEMS

A SISO fractional system is governed by a fractional differential equation:

y (t) + a 1 D α 1 y (t) + • • • + a m A D α m A y (t) = b 0 D β 0 u (t) + b 1 D β 1 u (t) + • • • + b m B D β m B u (t)
where (a j , b i ) ∈ R 2 , and the differentiation orders α 1 < α 2 < . . . < α m A , β 0 < β 1 < . . . < β m B are allowed to be non-integer positive numbers. State space representation was extended by [START_REF] Matignon | Some results on controllability and observability of finitedimensional fractional differential systems[END_REF] to commensurate fractional systems, where all the differentiation orders are multiple integers of α. The extension was done by allowing the differentiation order of the state vector to be any commensurate order α ∈ R + * . The fractional state space representation is presented in a MIMO case as:

D α x(t) = Ax(t) + Bu(t), (1)
y(t) = Cx(t) + Du(t) (2)
where x ∈ R n is the state vector, u ∈ R m the input vector, y ∈ R p the output vector,

A ∈ R n×n , B ∈ R n×m , C ∈ R p×n , D ∈ R p×m
are constant matrices. Here, zero initial conditions are considered: x(t) = 0 for t ≤ 0. Matignon [START_REF] Matignon | Stability properties for generalized fractional differential systems[END_REF] proved that the fractional system (1)-( 2) is stable if:

0 < α < 2 and | arg(λ k )| > α π 2 ∀k = 1, . . . , n
where λ k is the k th -eigenvalue of A and -π < arg(λ k ) ≤ π.

The conversion of ( 1)-( 2) to the MIMO transfer function form is obtained as for the rational systems by:

G(s) = C(s α I -A) -1 B + D
where s is the Laplace variable.

In the following, the pair (A, B) is assumed to be reachable and the pair (C, A) is assumed to be observable. The controllability and the observability conditions of a state space representation of a commensurate fractional system are the same as for rational systems [START_REF] Matignon | Some results on controllability and observability of finitedimensional fractional differential systems[END_REF].

One of the main difficulties with fractional models is the time-domain simulation. This problem has been extensively studied and an overview of the principal methods can be found in [START_REF] Aoun | Numerical simulations of fractional systems: an overview of existing methods and improvements[END_REF]. The most commonly used approximation of fractional operators is the recursive distribution of zeros and poles, proposed in [START_REF] Oustaloup | La dérivation non entière. Théorie, synthèse et applications[END_REF], which approximates the frequency behavior of s α in the frequency range [ω A , ω B ]. Nevertheless, this approximation has null asymptotic behaviors at low and high frequencies, which can introduce a static error between the fractional model and its approximation. To avoid this drawback, Trigeassou et al. [START_REF] Trigeassou | Modeling and identification of a non integer order system[END_REF] suggested to use the conventional integrator outside the frequency range [ω A , ω B ]:

s -α [ω A ,ω B ] = G α s N c ∏ k=1 1 + s/ω ′ k 1 + s/ω k (3)
where:

• N c is the number of cells (directly related to the quality of the approximation), • G α is fixed so that s -α has the same gain as

s -α [ω A ,ω B ] in the middle of the interval [ω A , ω B ],
• ω ′ k and ω k are respectively zeros and poles recursively distributed in the frequency range

[ω b , ω h ] = [σ -1 ω A , σω B ]
where σ is generally set to 10 to minimize border effects. They are defined by the following relations:

ω ′ k = γω k , ω k+1 = ηω ′ k , α = 1 - log γ log γη .
This approximation is used to simulate the fractional systems presented in this paper (section 5) with the parameters: N c = 20, ω A = 10 -5 and ω B = 10 5 .

SUBSPACE ALGORITHMS FOR FRACTIONAL STATE-SPACE IDENTIFICATION

Consider the linear continuous-time fractional state-space representation (1)-(2) corrupted by additive noise:

D α x(t) = Ax(t) + Bu(t) + v(t), (4)
y(t) = Cx(t) + Du(t) + w(t) (5)
where v ∈ R n and w ∈ R p are zero-mean processes. Under some mild conditions on v and v (uniform spectrum, uncorrelation) [START_REF] Johansson | Stochastic theory of continuous-time state-space identification[END_REF], the state-space representation (4)-( 5) can be replaced by the innovations model:

D α x(t) = Ax(t) + Bu(t) + Ke(t), (6)
y(t) = Cx(t) + Du(t) + e(t). (7)
The problem in this section is to estimate the system matrices A, B, C, D from sampled input-output data1 {u k } N-1 k=0 and {y k } N-1 k=0 . The commensurate order α is assumed to be known. The case where α is unknown is discussed in section 4.

To introduce the difficulty of the continuous-time fractional state-space model identification, let us consider firstly the deterministic case (e(t) = 0). Then, by computing the successive α-order fractional derivatives of [START_REF] Aoun | Synthesis of fractional Laguerre basis for system approximation[END_REF] and by substitution, the following extended linear model is obtained:

ȳ(t) = Γ * i x(t) + Φ * i ū(t) (8)
with input and output variables:

ū(t) = [u(t) T D α u(t) T . . . D iα u(t) T ] T ȳ(t) = [y(t) T D α y(t) T . . . D iα y(t) T ] T
and

Γ * i =      C CA . . . CA i-1      ∈ R ip×n , Φ * i =       D 0 • • • 0 CB D . . . . . . . . . . . . . . . 0 CA i-2 B • • • CB D       ∈ R ip×im .
The structure of ( 8) is the same than the extended linear model used in classical discrete-time subspace identification methods [START_REF] Katayama | Subspace methods for system identification[END_REF]. Unfortunately, (8) contains the successive α-order fractional derivatives of the input-ouput data which are not measured in most practical cases and which are difficult to estimate particularly in a noisy framework [START_REF]Identification of Continuous-Time Models from Sampled Data[END_REF].

To avoid this difficulty, the following operator (fractional low-pass filter) is introduced:

Λ(s) = 1 1 + s ω f α = 1 1 + τs α with τ = (1/ω f ) α . ( 9
)
Let us consider the Laplace transform of ( 6)-( 7):

s α X(s) = AX(s) + BU(s) + KE(s) ( 10
)
Y (s) = CX(s) + DU(s) + E(s). ( 11
)
Then, ( 10) can be expressed as:

X(s) = (I + τA)[Λ(s)X(s)] + τB[Λ(s)U(s)] + τK[Λ(s)E(s)] = A λ [Λ(s)X(s)] + B λ [Λ(s)U(s)] + K λ [Λ(s)E(s)]
with A λ = I + τA, B λ = τB and K λ = τK. Application of the inverse Laplace transform leads to the following system of linear equations:

x(t) = A λ [λx(t)] + B λ [λu(t)] + K λ [λe(t)] (12) y(t) = Cx(t) + Du(t) + e(t) (13)
where λx(t), λu(t) and λe(t) correspond to the states, the inputs and the noise prefiltered by Λ in [START_REF] Malti | Advances in system identification using fractional models[END_REF]. Then, from [START_REF] Johansson | Stochastic theory of continuous-time state-space identification[END_REF], it is found by recursion that:

y(t) = Cx(t) + Du(t) + e(t) = CA λ [λx(t)] + CB λ [λu(t)] + Du(t) + CK λ [λe(t)] + e(t) . . . = CA k λ [λ k x(t)] + k ∑ j=1 CA k-j λ B λ [λ k-j+1 u(t)] + Du(t) + k ∑ j=1 CA k-j λ K λ [λ k-j+1 e(t)] + e(t)
for k ∈ N * , where λ k x(t) denotes the signals obtained from x(t) by filtering through a series of k low-pass filters Λ. In the same way, it is found for l ∈ N * :

λ l y(t) = C[λ l x(t)] + D[λ l u(t)] + [λ l e(t)] = CA λ [λ l+1 x(t)] + CB λ [λ l+1 u(t)] + D[λ l u(t)] + CK λ [λ l+1 e(t)] + [λ l e(t)] . . . = CA k-l λ [λ k x(t)] + k-l ∑ j=1 CA k-j-l λ B λ [λ k-j+1 u(t)] + D[λ l u(t)] + k-l ∑ j=1 CA k-j-l λ K λ [λ k-j+1 e(t)] + [λ l e(t)]
with k ≥ l. As a consequence, the input-output data can be formulated as the following extended linear model (with no time derivatives of the data):

Y (t) = Γ i X (t) + Φ i U (t) + Ψ i E (t) (14)
with state, input-output and noise variables:

X (t) = λ i-1 x(t) Y (t) =       λ i-1 y(t) λ i-2 y(t) . . . λ 1 y(t) y(t)       , U (t) =       λ i-1 u(t) λ i-2 u(t) . . . λ 1 u(t) u(t)       , E (t) =       λ i-1 e(t) λ i-2 e(t) . . . λ 1 e(t) e(t)      
and:

Γ i =      C CA λ . . . CA i-1 λ      , Φ i =       D 0 • • • 0 CB λ D . . . . . . . . . . . . . . . 0 CA i-2 λ B λ • • • CB λ D       , Ψ i =       I 0 • • • 0 CK λ I . . . . . . . . . . . . . . . 0 CA i-2 λ K λ • • • CK λ I      
where Γ i ∈ R ip×n is the extended observability matrix and Φ i ∈ R ip×im , Ψ i ∈ R ip×im are block Toeplitz matrices. Now, from N available input-output samples observed at discrete times t k = kT s for k = 0, . . . , N -1, the extended linear model ( 14) can be rewritten as:

Y N = Γ i X N + Φ i U N + Ψ i E N ( 15
)
where

U N =        [λ i-1 u] 0 [λ i-1 u] 1 • • • [λ i-1 u] N-1 [λ i-2 u] 0 [λ i-2 u] 1 • • • [λ i-2 u] N-1 . . . . . . . . . [λu] 0 [λu] 1 • • • [λu] N-1 u 0 u 1 • • • u N-1        ∈ R mi×N and [λ j u] k = λ j u(t k ) denotes the sampled filtered data. The matrices Y N ∈ R pi×N , X N ∈ R n×N and E N ∈ R mi×N are con-
structed in a similar way. The formulation given in [START_REF] Ohsumi | Subspace identification for continuous-time stochastic systems via distribution-based approach[END_REF] enables to use subspace identification algorithms as in their original nonfactional discrete-time version. The difference is the addition of a step in which the data are filtered (what is a classical step in continuous-time identification [START_REF]Identification of Continuous-Time Models from Sampled Data[END_REF]). This lowpass filtering induces the tuning of an additional parameter ω f whose the sensibility has been studied in [START_REF] Thomassin | Subspace method for continuous-time fractional system identification[END_REF].

Deterministic framework : MOESP algorithm First, let us start by using the most popular subspace identification method called MOESP (MIMO output-error state space) algorithm [START_REF] Katayama | Subspace methods for system identification[END_REF][START_REF] Viberg | Subspace-based methods for the identification of linear time-invariant systems[END_REF]. This is a deterministic approach based on the properties of the noiseless version of (15):

Y N = Γ i X N + Φ i U N . ( 16
)
The principle of this algorithm is as follows:

1. Compute the LQ decomposition of the data matrix:

U N Y N = L 11 0 L 21 L 22 Q T 1 Q T 2 ( 17
)
where

L 11 ∈ R im×im , L 21 ∈ R ip×im , L 22 ∈ R ip×ip with L 11 ,
L 22 lower triangular, and

Q 1 ∈ R N×im , Q 2 ∈ R N×ip are or- thogonal. 2.
Compute the singular value decomposition (SVD) of the L 22 matrix approximating the column space of Γ i :

L 22 = U 1 U 2 Σ 1 0 0 0 V T 1 V T 2 ( 18
)
where

U 1 ∈ R ip×n , U 2 ∈ R ip×(ip-n) and Σ 1 ∈ R n×n
6. Estimate the B λ and D matrices. For that purpose, it can be shown that:

U T 2 Ψ i = U T 2 L 21 L -1 11 ( 20
)
which is a linear equation with respect to B λ and D. Define:

U T 2 L 1 L 2 • • • L i (21) U T 2 L 21 L -1 11 M 1 M 2 • • • M i (22) with L k ∈ R (ip-n)×p and M k ∈ R (ip-n)×m for k = 1, . . . , i.
Thus, from (20):

L 1 D + L 2 ĈD + • • • + L i Ĉ Âi-2 λ B λ = M 1 L 2 D + L 3 ĈD + • • • + L i Ĉ Âi-3 λ B λ = M 3 . . . L i-1 D + L i ĈD = M i-1 L i D = M i . Define L k = [L k . . . L i ] ∈ R (ip-n)×(i+1-k)p , k = 2, .
. . , i, and get the following overdetermined system of linear equations:

       L 1 L 2 Γi-1 L 2 L 3 Γi-2 . . . . . . L i-1 L i Γ1 L i 0        D B λ =        M 1 M 2 . . . M i-1 M i        (23)
where the block coefficient matrix in the left-hand side is

i(ip -n) × (p + n)-dimensional.
Estimates of B λ and D are obtained by finding the least-squares solution of [START_REF]Identification of Continuous-Time Models from Sampled Data[END_REF].

The matrices of the fractional continuous-time state-space representation ( 10)-( 11) are then deduced as follows: Â = 1 τ ( Âλ -I), B = 1 τ Bλ . The Ĉ and D matrices do not change. This method is consistent in its discrete-time version if the output noise (w) is white and if there is no process noise (v = 0). In our framework, the filtering of the input-output data calls the consistency into question.

Stochastic framework: PO-MOESP algorithm

In a stochastic framework, instrumental variables can be used to remove the noise effect. The most popular method in this context is the PO-MOESP (Past Output MIMO Output-Error State Space) algorithm [START_REF] Viberg | Subspace-based methods for the identification of linear time-invariant systems[END_REF][START_REF] Verhaegen | Identification of the deterministic part of MIMO state space models given in innovations form from input-output data[END_REF]. First, the input and output data and the noise matrices are partitioned into two parts:

U N = U p U f , Y N = Y p Y f , E N = E p E f . ( 24
)
The number of block rows in U p , Y p and E p is denoted by β. The number of block rows in U f , Y f and E f is then γ = i -β. These numbers can be chosen arbitrary, but it is generally preferred to set β = γ. By considering only the second part of these matrices, [START_REF] Ohsumi | Subspace identification for continuous-time stochastic systems via distribution-based approach[END_REF] becomes:

Y f = Γ γ X f + Φ γ U f + Ψ γ E f . ( 25
)
The instrumental variable introduced to remove the effects of the

noise is Z = [U T p Y T p ] T .
Hence, the principle of the PO-MOESP algorithm is as follows:

1. Compute the LQ decomposition of the data matrix:

    U f U p Y p Y f     =     L 11 0 0 0 L 21 L 22 0 0 L 31 L 32 L 33 0 L 41 L 42 L 43 L 44         Q T 1 Q T 2 Q T 3 Q T 4     (26)
where all the matrices L j j are lower triangular and all the matrices Q j are orthogonal (for j = 1, 2, 3, 4). 2. Compute the SVD of [L 42 L 43 ] which approximates the column space of Γ γ :

L 42 L 43 = U 1 U 2 Σ 1 0 0 0 V T 1 V T 2 ( 27
)
where

U 1 ∈ R ip×n , U 2 ∈ R ip×(ip-n) and Σ 1 ∈ R n×n .
As for the MOESP algorithm, the state-space order n can be estimated by calculating the number of significant singular values of L 42 L 43 .

The following steps 3, 4, 5 and 6 are the same as for the MOESP algorithm.

FRACTIONAL ORDER ESTIMATION

In this section, we assume that the fractional derivative order α ∈ (0, 2) is unknown and has to be estimated by minimizing a

1 2 ŷc (α) -y c 2 2 , (28)
where y c is the vector (of length pN) obtained by concatenating the p system outputs and ŷc (α) is the vector obtained by concatenating the p outputs of the state-space representation estimated with one of the subspace methods (MOESP or PO-MOESP) for a given α. The value of α in ( 28) is obtained by a nonlinear optimization technique. The iterative algorithm used is based on a subspace trust-region method and on the interior-reflective Newton method described in [START_REF] Coleman | An interior, trust region approach for nonlinear minimization subject to bounds[END_REF]. Consequently, the proposed subspace algorithm is executed at each iteration of the optimization technique. The number of state variables n is assumed to be known. However, it may be estimated by using order estimation criteria presented in [START_REF] Bauer | Order estimation for subspace methods[END_REF] .

SIMULATION EXAMPLES

The algorithms are applied to input-output data of length N = 1023 generated by simulating the linear system (1)-( 2) with α = 0.9,

A = 0 -0.1 1 -0.2 , B = 1 0 ,C = 0 0.1 0.5 -0.1 , D = 0 0 ,
and zero initial conditions (x(t) = 0 for t ≤ 0). The input signal is a pseudo-random binary sequence (PRBS) with maximum length. The sampling period is T s = 0.05s. The outputs (Fig. 1) are corrupted by white noise with a signal-to-noise ratio (SNR) of 20 dB. The influence of the method parameters (the number of block rows i and the filter frequency ω f ) has been studied in [START_REF] Thomassin | Subspace method for continuous-time fractional system identification[END_REF]. It was shown that ω f can be chosen in a suitable range over one decade, between 1 and 10 and that the choice of i influences more the results: acceptable results are obtained starting from i = 6. For this example, the algorithm parameters are fixed to i = 8 and ω f = 6.

Using the MOESP algorithm

The estimation results with the MOESP algorithm are: To analyze the statistical properties of the estimator, three Monte Carlo simulations are carried out for SNR of 20, 15 and 10 dB. Each simulation is done with 500 runs with different realizations of noise. The means of the normalized error norms, obtained with i = 8 and ω f = 6, are given in Tab. 1. Figure 2 shows the estimated poles. The normalized mean squared error (MSE) of the poles are indicated in Tab. 2. It can be seen that the pole estimator is biased. Moreover, as shown in the histograms of the estimated order α in Fig. 3, the order estimate is also biased, and the bias increases with the noise level. These biases are obviously introduced by the MOESP algorithm. Indeed, in order to guarantee the consistency of the estimates using this algorithm, the input must satisfy the persistent excitation condition (verified with a PRBS for example) and the output noise have to be white [START_REF] Bauer | Analysis of the asymptotic properties of the MOESP type of subspace algorithms[END_REF]. In fact, the latter condition is not satisfied because of the data filtering.

α = 0.8894, Â = 0 -0.1008 1 -0.1757 , B = 1 0 , Ĉ = 0.

Using the PO-MOESP algorithm

The PO-MOESP algorithm applied to the same input-output data gives the following estimates: The normalized prediction error norm equals 0.0983 (-20.15 dB) for the first output and 0.0989 (-20.09 dB) for the second one. So the results are improved with the PO-MOESP algorithm. The same Monte Carlo simulations are performed as previously. Table 1 gives the means of the normalized error norms. Figure 4 shows the estimated poles. Their normalized MSE are indicated in Tab. 2. The estimates are unbiased but a slight increase of the variance is observed as compared to the results obtained with the MOESP algorithm. However, the normalized MSE is smaller, especially for a low SNR.

As shown in the histograms of the estimated order α in Fig. 5, the estimation is now unbiased. Moreover, the variance on the estimated parameter α is improved.

CONCLUSION

This paper focuses on the identification of continuous-time fractional state-space representation. Thanks to an adapted data filtering, we have shown that subspace-based algorithms can be used to estimate the state-space matrices. Results on the identification of a multivariable system with the well-known MOESP and PO-MOESP algorithms are given. The commensurate dif- ferentiation order is estimated by using nonlinear programming. Simulation examples have shown that the estimators are biased with the MOESP algorithm. However, the PO-MOESP algorithm eliminates this bias by using an instrumental variable. The limitation of this approach is that the design of the instruments requires an increased number of block rows i, i.e. an increase of the series of fractional low-pass filters Λ(s), which leads to an excessive computational cost. So, a future work could be to consider another instruments to solve this problem.

4 .

 4 . The state order n can be estimated from the SVD since n = dim Σ 1 in the noiseless case. 3. Estimate the extended observability matrix: Γi = U 1 Σ Estimate the C matrix: Ĉ = Γi (1 : p, 1 : n).

5 .

 5 Estimate A λ by solving the linear equation: Γi 1 : p(i -1), 1 : n A λ = Γi p + 1 : ip, 1 : n .

Figure 1 .

 1 Figure 1. OUTPUT DATA y AND OUTPUT OF THE ESTIMATED MODEL ŷ WITH MOESP (SNR=20 dB, i = 8 AND ω f = 6).

  0.1007 0.5008 -0.1020 , D = 0.4564 10 -3 0.0491 10 -3 .

Figure 2 .

 2 Figure 2. ESTIMATED POLES WITH MOESP, WHERE + DENOTES THE TRUE POLES.

Figure 4 .

 4 Figure 4. ESTIMATED POLES WITH PO-MOESP, WHERE + DE-NOTES THE TRUE POLES.

Figure 5 .

 5 Figure 5. HISTOGRAMS OF THE ESTIMATED COMMENSURATE ORDER α WITH PO-MOESP.

  The normalized prediction error norm ŷy 2 / y 2 equals 0.1050 (-19.57 dB) for the first output and 0.1008 (-19.93 dB) for the second one.

	0012 0.1012 0.4997 -0.0923	, D =	1.3782 10 -3 -0.6763 10 -3 .

Table 1 .

 1 MEANS OF THE NORMALIZED ERROR NORMS.

		output SNR=20 dB SNR=15 dB SNR=10 dB
	MOESP	1 2	0.1026 (-19.78 dB) (-14.73 dB) (-9.56 dB) 0.1834 0.3327 0.1005 0.1797 0.3282
			(-19.95 dB) (-14.91 dB) (-9.68 dB)
	PO-MOESP	1 2	0.1000 (-20.00 dB) (-15.09 dB) (-10.35 dB) 0.1760 0.3038 0.0991 0.1745 0.3010 (-20.08 dB) (-15.16 dB) (-10.43 dB)

Table 2 .

 2 NORMALIZED MSE OF THE POLES.

SNR

MOESP PO-MOESP 20 dB (-0.4952 ± i 0.8875).10 -4 (0.4661 ± i 0.5006).10 -4

15 dB (-0.2530 ± i 0.6315).10 -3 (0.1634 ± i 0.1976).10 -3 10 dB (-0.1208 ± i 0.5740).10 -2 (0.4651 ± i 0.7347).10 -3

  Figure 3. HISTOGRAMS OF THE ESTIMATED COMMENSURATE ORDER α WITH MOESP.

		60									60									60
		50									50									50
	rate (in %)	30 40								rate (in %)	30 40								rate (in %)	30 40
		20									20									20
		10									10									10
		0.82 0	0.84	0.86	0.88	0.9	0.92	0.94	0.96	0.82 0	0.84	0.86	0.88	0.9	0.92	0.94	0.96	0.82 0	0.84	0.86	0.88	0.9	0.92	0.94	0.96
				commensurate order	α					commensurate order	α					commensurate order	α
				(a) SNR=20 dB							(b) SNR=15 dB							(c) SNR=10 dB
		60									60									60
		50									50									50
	rate (in %)	30 40								rate (in %)	30 40								rate (in %)	30 40
		20									20									20
		10									10									10
		0.82 0	0.84	0.86	0.88	0.9	0.92	0.94	0.96	0.82 0	0.84	0.86	0.88	0.9	0.92	0.94	0.96	0.82 0	0.84	0.86	0.88	0.9	0.92	0.94	0.96
				commensurate order	α					commensurate order	α					commensurate order	α
				(a) SNR=20 dB							(b) SNR=15 dB							(c) SNR=10 dB

The discrete-time variables are denoted by x k and correspond to the time sampling with a constant sampling period T s of the continuous-time variable x(t): x k = x(kT s ).

Copyright c 2009 by ASME


Jérôme Lacan

Tanguy Pérennou

Reducing Retransmissions in Point-to-point Satellite Transmissions

This paper presents a mechanism that reduces the number of retransmissions needed to ensure reliable unicast transmissions over a shared broadcast channel. This improves the bandwidth occupancy by information packets and potentially reduces the end-to-end jitter and delay. After each block transmission, the mechanism determines the minimal set of information packets to retransmit, builds a specific Reed-Solomon code from this set and then sends repair packets produced with this code instead of retransmitting the information packets. The mechanism was analysed theoretically and evaluated numerically under many different conditions (variable packet loss rate, number of receivers, etc.). Both the erasure channel and the error channel were explored. In most cases, the system can be set up such that it uses less retransmissions than a classical ARQ mechanism: most of the time, 40% up to 95% less. Such a mechanism is clearly useful in satellite contexts.

I. Introduction

This paper deals with reliable unicast transmissions using a satellite communication system. Satellite communication systems mainly use the radio channel to send information packets to terminals or to exchange other control or management packets. As many other wireless systems (2.4GHz and 5GHz Wireless LANs, WiMax, cellular phone systems, etc.), this broadcast channel can have a high bit error rate (BER) that leads to a lot of packet losses as compared to wired communication systems. In addition, they exhibit a comparatively low data rate. Improving the goodput of such systems is therefore a useful task.

A first solution consists in reducing the BER of the channel using improved antennas and signal processing. Additionally, the reliability of satellite transmissions is classically ensured by pro-active mechanisms such as error correcting codes at the physical layer. Indeed, contrary to Wireless LANs, reactive mechanisms such as the Automatic Repeat reQuest (ARQ) are generally not implemented at the link layer because the large round-trip time on the satellite link can lead to unacceptable delays (especially with a geostationary satellite) when several consecutive retransmissions are needed.

If ARQ does not seem appropriate to improve the satellite network performance, error correcting codes or erasure codes (also called Forward Error Correction (FEC)) provide a well-known alternative to retransmissions by generating a set of repair packets from a set of information packets. Both information and repair packets are sent on the channel, and if some information packets are lost or received with bit errors, they can be recovered using the received repair packets. In some cases, the joint use of retransmission techniques and error correcting techniques, called Hybrid-ARQ (H-ARQ), allows to combine the advantages of these two techniques (see e.g. Ref. 1).

In this paper we propose to use a generalization of H-ARQ techniques at the data link level, in such a way that each packet retransmitted on the channel can recover other lost or erroneous packets in several concurrent point-to-point transmissions. We first consider a packet erasure channel in which all packets containing errors are erased and do not reach the receiver: such packets lost in concurrent transmissions are considered as information packets in an erasure code. Each retransmitted packet is then a repair packet generated from the lost information packets. We then generalize the proposed mechanism to an error channel in which packets containing errors do reach the receiver and are considered as information packets in an error correcting code. Similarly, each retransmitted packet is a repair packet generated from the erroneous information packets. The paper will demonstrate theoretically and numerically that the amount of retransmissions dramatically decreases as compared to the classical ARQ approach.

The paper first reviews related work (Section II). The principle of our proposition is then detailed (Section III). Potential gains are then evaluated in a theoretical manner (Section IV) and numerical results are presented and discussed (Section V). Then Section VI shows how the proposed mechanisms can be applied to several satellite communication contexts. Finally, Section VII concludes and presents future work in this direction.

II. Related Work

In the general context of this paper, several users share a broadcast channel for point-to-point transmissions. We assume that the channel can only broadcast a packet at a time and that transmissions are temporally multiplexed at the packet level. This type of channel is used in wireless local or metropolitan area networks (Bluetooth, 802.11, WiMAX), mobile networks (GSM, GPRS, UMTS) or satellite networks.

The above types of network rely on the radio resource, which is particularly characterized by a low transmission rate and a high bit error rate (BER) as compared to wired networks. Error control mechanisms are used at the physical level and allow the correction of a part of bit errors.

Other mechanisms allow the correction of the remaining errors so as to ensure end-to-end reliability. In the case of the Transmission Control Protocol (TCP) 2 for instance, information packets that are lost are not acknowledged by the receiver and several successive retransmissions are attempted upon expiration of a timer. As TCP was designed to handle congestions for wired networks where router congestions are the principal cause of packet losses, it is not optimal when packet losses are due to bit errors. TCP's additive increase multiplicative decrease (AIMD) congestion control reacts to a packet loss by both a decrease of the sending rate and a doubling of the retransmission time-out. If the loss is due to a bit error, those actions are useless and simply decrease bandwidth occupancy although there is no congestion in the network. [START_REF] Caceres | Improving the Performance of Reliable Transport Protocols in Mobile Computing Environments[END_REF][START_REF] Balakrishnan | A Comparison of Mechanisms for Improving TCP Performance over Wireless Links[END_REF][START_REF] Holland | Analysis of TCP Performance Over Mobile Ad Hoc Networks[END_REF] The same problem exists for more other transport protocols using congestion control like TFRC, 6 DCCP 7 or SCTP. [START_REF] Stewart | Stream Control Transmission Protocol[END_REF] Several propositions have been made to implement end-to-end loss differentiation algorithms and improve end-to-end performance. [START_REF] Cen | End-to-end differentiation of congestion and wireless losses[END_REF] In general wireless networks, retransmission mechanisms are used at the data link layer to recover as much as possible packet losses due to bit errors. In the 802.11 family MAC layer 10 for example, each data frame may be retransmitted several times until it is acknowledged by the receiver; if several retransmissions in a row fail, the auto-rate fallback mechanism is activated and switches to a modulation that is more robust to error bits (at the expense of a lower transmission rate). Retransmissions are also used in the Radio Link Control Protocol (RLC) layer of UMTS in its acknowledgement mode [START_REF] Gpp | Radio Link Control (RLC) Protocol Specification[END_REF] as well as in recent propositions for satellite link layers. [START_REF] Kojo | Improving TCP Performance over Wireless WANs using TCP/IP-Friendly Link Layer[END_REF] The use of retransmissions at the lower layers however causes a number of problems. When the number of retransmissions gets too high, the useful bandwidth occupancy decreases while end-to-end delays increase. When those delays increase too much, higher level error recovery mechanisms such as those cited above for TCP are automatically unleashed, and the overall network performance unavoidably decreases. [START_REF] Fairhurst | Advice to Link Designers on Link Automatic Repeat reQuest (ARQ)[END_REF] On the erasure channel, a well-known alternate solution is the use of an erasure code which allows the correction of a certain amount of losses in a set of information packets by using repair packets. For instance, in the Maximum Distance Separation (MDS) family of erasure codes (including Reed-Solomon codes), the MDS(k, n) allows the computation of n -k repair packets for a block of k initial information packets. As soon as the receiver receives any k packets among the k +(n-k) packets, it will be able to decode the original k information packets. [START_REF] Rizzo | Effective Erasure Codes For Reliable Computer Communication Protocols[END_REF] Other families of erasure codes were defined. The most famous are Low-Density Parity-Check (LPDC) codes and Raptor codes which have very fast coding and decoding algorithms but which are not optimal in term of correction capability.

Most of these erasure codes (except Raptor codes) are derived from classical error correcting codes. The most famous are turbo codes, LDPC codes or algebraic codes like Reed-Solomon or BCH codes.

Erasure and error correcting codes can be used in a proactive way (Forward Error Correction or FEC) or in a reactive way (Hybrid ARQ). In the reactive way, the sender generates and sends redundancy following the feedbacks of the receiver.

Error correcting codes are generally used in proactive way, however, some of the most recent wireless standards such 3GPP HSDPA [START_REF] Gpp | Technical Specification Group Radio Access Network; High Speed Downlink Packet Access (HSDPA)[END_REF] or "Mobile Wimax" 802.16e [START_REF]IEEE Standard for Local and metropolitan area networks Part 16: Air Interface for Fixed and Mobile Broadband Wireless Access Systems Amendment 2: Physical and Medium Access Control Layers for Combined Fixed and Mobile Operation in Licensed Bands and Corrigendum 1[END_REF] define a framework to implement hybrid ARQ schemes. Erasure codes are mainly used for real-time transmissions and multicast/broadcast transmissions (see e.g. Ref. 17). In both applications, when a return channel is available, the level of redundancy can be adapted to te channel conditions so as to improve the global performance.

Recently, the concept of opportunistic coding was introduced to improve practical network coding performance. [START_REF] Katti | XORs in the air: practical wireless network coding[END_REF] The idea is to use some informations on the data needed by the neighbors (i.e the receivers) to determine how to combine the information packets in order to generate "optimal" redundant packets. In Ref. 19, this concept was applied to packet losses recovery in the context of unicast and multicast transmissions on a broadcast channel. The main result is the construction of an "optimal" repair packet with XOR operations.

Our proposition can be considered as a mix between opportunistic and hybrid ARQ techniques. To recover losses or errors of various unicast transmission on a broadcast channel, we first opportunistically determine the minimal set of needed information packets. Then, we generate the repair packets with classical but efficient erasure/error correcting codes. This proposal allows a dramatical reduction of the number of retransmissions compared to classical ARQ mechanism.

III. Applying Forward Error Correction Codes to Retransmissions

A. Principle

The mechanism proposed in this paper essentially aims at reducing the number of retransmissions at the data link level. This will allow both an increase of the proportion of information packets broadcasted on the channel, thus improving bandwidth occupancy, and as a side effect a decrease of the end-to-end delay and jitter.

We make the basic assumption that transmissions are all issued from a unique sender. This assumption corresponds for instance to satellite transmissions or transmissions relayed by an access point in a cellular network. We also make the assumption that the access to the broadcast channel is shared. The kind of multiplexing is generally dynamic, as in the Carrier Sense Multiple Access (CSMA) found in Ethernet and 802.11, or in the Demand-Assigned Multiple Access (DAMA) used in DVB-RCS return link of satellite communications.

We adapt the Hybrid ARQ type II mechanism 1 to point-to-point transmissions by coupling the filtering function and the retransmission mechanism at the data link level.

The proposed mechanism is based on the fact that on a broadcast channel each receiver receives at the physical level all the broadcasted packets (promiscuous reception). In our proposition, instead of filtering out packets according to their destination address, each receiver temporarily archives all information packets received at the physical level. Then each receiver acknowledges only the packets addressed to him. Then the sender computes the list of non acknowledged packets and builds repair packets based on that list. The list of unacknowledged packet IDs is broadcasted to all receivers, then the sender broadcasts redundancy packets until all receivers acknowledge all the previously lost or rejected information packets. A receiver can acknowledge such a packet as soon as he has received enough redundancy: in the case of a packet erasure channel, he must have received at least as many repair packets as lost packets in the list, and for an error channel, he must have received enough repair symbols to correct erroneous information symbols in the list.

B. Illustration

Figure 1 illustrates the proposition for an erasure channel. First, sender S, typically a satellite, relays information packets to receivers R 1 and R 2 using a simple interleaving scheme. Each receiver should receive four information packets for himself. We note info(i, j) the i th information packet sent to the j th receiver. R 1 does not receive his first information packet while R 2 does not receive his second and third information packets. Both receivers acknowledge all correcly received packets that were addressed to them to S. S can then establish that the list of unacknowledged information packets IDs is L = {(1, 1); (2, 2); (3, 2)}. S broadcasts the L list to R 1 and R 2 , then builds a Reed-Solomon code based on the information packets of this list (Reed-Solomon codes are not the only option, especially when this mechanism is adapted to an error correcting channel).

S then sends repair packets made with the previously built code until all receivers have acknowledged their four packets. A first repair packet -repair(1) -is broadcasted to both receivers, which receive it without error. Now R 1 has buffered packets info(2,2), info (3,2) and repair(1), so he can decode the L list and recover missing packet info (1,1). R 1 then acknowledges its block to S. R 2 , however, misses information packets info(2,2) and info(3,2) from L and even with repair(1) he still does not have enough redundancy to decode L. S has not received the block acknowledgement from R 2 and therefore sends repair (2), which is lost by R 2 . S still does not receive any block acknowledgement from R 2 and sends repair(3) which is correctly received by R 2 . With info(1,1), repair(1) and repair(3), R 2 decodes the L list and recovers missing packets info(2,2) and info (3,2), and acknowledges its block to S.

In this example, only 3 retransmissions were needed for both users to get their information packets. With a classical ARQ mechanism, 4 retransmissions would have been needed (the second retransmitted packet is lost and needs additional retransmission).

C. Sender and Receiver Algorithms

We consider an erasure channel and we assume a static temporal multiplexing: if R receivers simultaneously receive their own N packets at the data link level, and if we note info(i, j) the transmission of the i th information packet addressed to the j th receiver on the broadcast channel, then the order of the interleaved information packets on the channel will be info(1,1), ..., info(1,R), info(2,1), ..., info(N ,R).

The sender's algorithm is basically as follows:

1. send a block made of N × R interleaved information packets (B n );

2. determine the complete list of packets lost by their addressee in this block (L n );

3. consider the packets in L n as the k information packets of a Reed-Solomon (RS) code and build code RS n ;

4. broadcast the list of the L n packet IDs;

5. broadcast repair packets produced by the RS n code until all receivers have sent an acknowledgement for the full block B n .

Concerning receivers, each one promiscuously listens to the broadcast channel and uses the following algorithm:

1. temporarily archive all packets received in block B n including the packets not explicitly addressed to him;

2. acknowledge each received packet that is addressed to him;

3. receive the list of packet IDs forming the new RS n code;

4. receive and archive the successive repair packets of the RS n code;

5. as soon as there are enough repair packets to decode L n , recover lost packets addressed to him from the decoded L n and acknowledge block B n .

For the error channel, the sender's algorithm is the same. On the receiver side, each receiver acknowledges only the error-free packets addressed to him. He stores all received packets, erroneous or not, so as to exploit as much valid symbols as possible. Upon reception of each repair packet, the receiver tries to decode the erroneous packets of the list L n with the set of received repair packets. When the decoding is successful, it sends a block acknowledgment to the sender.

In both erasure and error channels, the repair packets are built from a Reed-Solomon (see Figure 2). For the erasure channel, the packets are built following the classical construction 14 by considering the packets as columns of a matrix and by encoding the rows. For the error channel, each packet is virtually split into G groups of b symbols. The m repair packets are also made of G groups of repair symbols, such that each "row" made of (N.R + m) × b symbols is a code-word with N.R.b information symbols and mb repair symbols. Each row can be then decoded independently from the others.

IV. Analytical Evaluation

In this Section we provide an analytical evaluation of the mechanism described in Section III in the case of both the erasure channel and the error channel. We use the following assumptions:

• Homogeneous receivers: all R receivers have the same properties.

• Same number of information packets per receiver per block: each receiver is addressed N packets in each block.

• Fixed-sized packets: all packets have the same size. They are made of G groups of b symbols. These groups may be used to perform independent error correction attempts on different groups.

• Independent erroneous packets: the distribution of erroneous packet does not depend on time nor on receivers (no temporal nor spatial correlation). The probability that a packet contains errors is p and therefore the probability that it does not contain any error is 1 -p. On an erasure channel, erroneous packets are considered lost. On an error channel, erroneous packets are considered received, but with a subset of erroneous symbols.

• Independent symbol errors: when a packet is erroneous, its errors are uniformly distributed with a probability q (symbol error rate) over the G.b symbols of the packet.

The remainder of this Section evaluates the overhead of the proposed mechanism for the erasure channel and the error channel, measured as the amount of retransmissions needed for each receiver to recover its N information packets.

A. ARQ Overhead

In a classical ARQ mechanism 1 each erroneous packet is retransmitted until reception without error. The overhead is simply computed as the mean number of retransmissions for one packet multiplied by the number of packets. This overhead is denoted O ARQ .

Let A be the random variable (r.v.) representing the number of trials needed to successfully transmit a single packet. Under the above uniform independent assumptions, we have P (A = a) = p (a-1) (1 -p) and the amount of trials needed is E[A] = 1/(1 -p). The amount of retransmissions needed for one packet is 1 -E[A], and therefore:

O ARQ = N.R.p 1 -p (1)
This value will be used as a basis of comparison for the numerical evaluation of our mechanism in Section V.

B. Hybrid ARQ Proposal Overhead

Let M be the r.v. representing the minimal amount of repair packets such that all R receivers can decode the N information packets addressed to them. The overhead of the proposed hybrid ARQ mechanism is therefore the expectation of M :

O HARQ = E[M ].
Let L be the r.v. representing the size of the list of erroneous packets discussed in Section III. To compute P (M = m), we successively switch to the CDF of M , introduce the possible values of L and apply Bayes' theorem:

P (M ≤ m) = R.N =0 P (M ≤ m | L = ).P (L = ) (2)
L follows a binomial law B(N.R, p) :

P (L = ) = N.R p (1 -p) N.R- (3)
Let S be the r.v. representing the amount of repair packets needed by a single receiver to decode all its information packets knowing that L = . The law of M is the maximum of the S j , j = 1, . . . , R of each receiver, which translates to the product of CDFs:

P (M ≤ m) = R.N =0 P (S ≤ m | L = ) R .P (L = ) (4)
To determine P (S ≤ m | L = ) for a given receiver, we introduce the following random variables: let X be the r.v. describing the amount of erroneous information packets that were addressed to him, let Y be the r.v. describing the amount of packets received with errors by this receiver in the L list (including the X erroneous packets), and let Z be r.v. describing the amount of repair packets received with errors by this receiver. Using Bayes again and exploiting independence properties, we can show that:

P (S ≤ m | L = ) = x,y,z P (S ≤ m | Y = y ∧ Z = z).P (X = x | L = ).P (Y = y | X = x ∧ L = ).P (Z = z) (5)
In the above equation, r.v. Z follows a binomial law B(m, p); r.v. Y with X = x and L = fixed can be determined knowing that Y -x follows a binomial law B( -x, p); and r.v. X with L = fixed follows a hypergeometric law H( , N, N.R -N ). The three last terms of Eq. 5 are therefore:

P (X = x | L = ) = N x . N.R -N -x / N.R ∀x, 0 ≤ x ≤ min( , N ) (6)
P (Y = y | X = x ∧ L = ) = -x y -x p y-x (1 -p) -y ∀y, x ≤ y ≤ (7)
P (Z = z) = m z p z (1 -p) m-z ∀z, 0 ≤ z ≤ m (8)
All other values of x, y and z lead to null terms.

Erasure Channel

With a Reed-Solomon code on the erasure channel, the receiver can recover erroneous packets erased by the channel if and only if he has received at least as many repair packets as erroneous packets (information or repair), i.e. m ≥ y + z. So when y and z are fixed:

P (S ≤ m | Y = y ∧ Z = z) = 1, if m ≥ y + z; = 0, otherwise. (9)
Hence for the erasure channel:

P (S ≤ m | L = ) = x,y,z P (X = x | L = ).P (Y = y | X = x ∧ L = ).P (Z = z) with: 0 ≤ x ≤ min( , N ), x ≤ y ≤ , 0 ≤ z ≤ max(0, m -y). (10)

Error Channel

On the error channel, we must consider the symbols in every packet, including erroneous information and repair packets. In Eq. 5, P (S ≤ m | Y = y ∧ Z = z) must be determined in the context of a Reed-Solomon error correcting code. The decoding is applied independently on the G code words constituted by the concatenation of all groups of b symbols located at the same offset in all stored information and repair packets. So here the information consists of groups of b symbols and the redundancy consists of m groups of b symbols. With a Reed-Solomon code, a code word can be decoded if and only if the amount of erroneous symbols is lesser or equal to half the amount of redundancy symbols, so at most m.b/2 erroneous symbols can be corrected. Let D be the r.v. such that P (D ≤ m) is the probability for a receiver to decode a code word with m or less repair packets. Under the independence assumptions above, the G code words can be decoded independently from each other, so:

P (S ≤ m | Y = y ∧ Z = z) = P (D ≤ m | Y = y ∧ Z = z) G (11)
Let E be the r.v. describing the total number of erroneous symbols in a single code word. As stated above, the code word can be decoded if and only if E ≤ m.b/2:

P (D ≤ m | Y = y ∧ Z = z) = P (E ≤ m.b/2 | Y = y ∧ Z = z) (12) = m.b/2 e=0 P (E = e | Y = y ∧ Z = z) (13)
Each code word contains (y + z)b symbols being part of erroneous packets. As symbol errors are independent in each erroneous packet, E follows a binomial law B((y + z)b, q):

P (E = e | Y = y ∧ Z = z) = (y + z)b e q e (1 -q) (y+z)b-e (14)
Hence for the error channel:

P (S ≤ m | L = ) = x,y,z   m.b/2 e=0 P (E = e | Y = y ∧ Z = z)   G × P (X = x | L = ).P (Y = y | X = x ∧ L = ).P (Z = z) with: 0 ≤ x ≤ min( , N ), x ≤ y ≤ , 0 ≤ z ≤ m. (15)

V. Numerical Evaluation

In this Section, we compare the overhead of the mechanism described in Section III with a classical ARQ retransmission mechanism. For both the erasure channel and the error channel, the evaluation is carried out using the equations given in Section IV. We define the gain of our mechanism over ARQ as the percentage of retransmissions needed by ARQ and not needed by HARQ:

gain = 100 × O ARQ -O HARQ O ARQ
A positive gain means that the proposed mechanism is better than classical retransmissions. An increasing gain denotes an increasing performance of HARQ, i.e. proportionally less and less retransmissions. The theoretical maximum gain of 100% is obtained when O HARQ = 0 and O ARQ > 0. A negative gain means that classical retransmissions are less numerous and behave better than our mechanism.

An extensive campaign of evaluation was made for both channels, varying all parameters: the packet error rate p, the size of symbol groups in a packet G such that G.b = 256 symbols, the number of receivers R and the amount of packets addressed to each receiver N . The symbol error rate q = 0.10 was fixed.

The best results are systematically obtained for G = 1 and b = 256. In that configuration, the error channel always allows a better gain than the erasure channel. Figure 3 shows the results for p = 0.01 (top figures) and p = 0.15 (bottom figures), N varying from 1 to 8 packets per block and R varying from 1 to 64 receivers. Both gains increase with N as well as with R. The error channel reaches a gain of 81% for p = 0.01 and of 95% for p = 0.15, while the erasure channel reaches a gain of 73% for p = 0.01 and of 72% for p = 0.15. These maximum values are all reached for the maximum N.R value. The Figure also shows that in most cases better gains are obtained for p = 0.15. In the case of the error channel, this is intuitively explained by the fact that there is a lot of erroneous packets but not a lot of erroneous symbols: a single packet containing 256 repair symbols can repair all the erroneous symbols needed by all the receivers.

The error channel does not always outperforms the erasure channel. The gain obtained with the error channel decreases as the number of groups per packet G increases from 1 to 256, and for many (N, R) couples there is a minimum value of G such that gain errors < gain erasure . This is particularily true for p = 0.01 where such a minimum G exists for all investigated (N, R) couples, as illustrated in Figure 4. The minimum G value then varies between 32 and 128 (we did not perform an exhaustive campaign of measure and sticked to powers of 2 for N , R and G, so the actual threshold values may be lower than that). Going beyond that G value leads to a dramatic decrease of the error channel performance. In the worst case, i.e. N = 1 and R = 2, setting G = 256 on the error channel leads our mechanism to use nearly 5 times more retransmissions than on the erasure channel.

To sum up, both proposed hybrid ARQ type II mechanisms provide potentially high gains on the number of retransmissions. The largest gains obtained correspond to higher values of the packet error rate (i.e. 15% rather than 1%) and large values of N and R. When a small amount of group symbols per packet can be chosen, the performance of the error channel is clearly better than the performance of the erasure channel.

Discussion

A number of constraints related to the mechanism proposed in this paper should be carefully examined.

One of these constraints is the storage capacity needed to temporarily store the packets addressed to other receivers. The needed storage capacity is a function of parameters R (number of receivers) and N (number of information packets per receiver and per block), as well as a function of the transmission rate. Yet, the higher N × R, i.e. the number of packets per block, the higher the efficiency of the proposed mechanism will be. Among these three parameters, only N and R may be adjusted by the sender so as to find a good compromise between the obtained gains (see Section V) and the amount of storage capacity. The number of receivers R may be controlled through the constitution of subgroups of jointly managed receivers. For instance on a broadcast channel shared among 10 receivers, 2 subgroups of 5 receivers can be formed, each receiver of a group ignoring the packets sent in the other subgroup. More sophisticated solutions may be elaborated, e.g. grouping receivers according to their storage capacity.

The number of information packets per receiver and per block N can be adjusted too. However augmenting N leads to an additional delay: in case of packet loss, the N information packets addressed to a given receiver will be delivered to the upper communication layers only after decoding the block. This can generate delay and jitter on the first information packets and make the use of e.g. the TCP protocol at the transport layer rather erratic, because the resulting delays may trigger unwanted retransmissions at the TCP level as well as a decrease of the transmission window [START_REF] Fairhurst | Advice to Link Designers on Link Automatic Repeat reQuest (ARQ)[END_REF] although the lost information packet will eventually be decoded by the data link layer. Yet, considering all N packets addressed to a receiver within a block, the last packet will be delivered faster than without the HARQ mechanism presented here. In practice the N parameter should be adjusted according to the maximum delay accepted by upper layers (especially when using TCP) and the storage capacities of the receivers.

Another potential implementation issue is the size of packets. A number of data link level protocols use a variable packet size, like the 802.11 protocol. The implementation should provide an operation that combines two packets of different sizes. In some very specific cases, for instance if a large packet is combined with several small packets, several large repair packets may be retransmitted and thus increase channel occupancy overhead. However this phenomenon can be partially avoided by the constitution of blocks with packets of the same size. The evaluation of the previous section assumes that packets losses are uniformly distributed over time, and that they are independent for different receivers. This assumption does not cover all kinds of broadcast channels. In particular in 802.11 wireless networks, packet losses are correlated in time and space, and it has been shown that these correlations have a strong impact on the performance of reliable multicast transmissions. [START_REF] Lacan | Evaluation of Error Control Mechanisms for 802.11b Multicast Transmissions[END_REF] A classical solution to cope with bursts of losses or errors is to use interleaving techniques. The mechanism proposed in this paper may therefore be evaluated under stronger assumption of spatial and temporal correlation in future work.

For the error channel, the choice of the parameter G is important. As shown on Figure 4, if the code has a limited length (for example, classical octet-based Reed-Solomon implementations do not work for a number of encoding symbol greater than 256), the parameter G should be set as small as possible. If, for the minimal value G, the decoding performances are not acceptable, other error correcting codes supporting a large number of encoding symbols (such e.g. LDPC codes) should be chosen.

VI. Applications in Satellite contexts

Several architectures can be considered for satellite communications. In this Section, we consider 3 classes of architecture: bi-directional point-to-point satellite communications, satellite broadcast with terrestrial repeaters, and satellite broadcast without repeaters, but with an ad-hoc network of receivers. For each architecture, we describe the general principle allowing to apply our proposal.

A. Satellite-only

In the satellite-only case, we consider that several bidirectional point-to-point communications share the same satellite link. A classical example is an IP communication using DVB-S or DVB-S2 on the forward link and DVB-RCS or a terrestrial link on the return link. The proposed mechanism can be implemented at the link layer between the gateway and the satellite terminals. It optimizes the link utilization by reducing significantly the number of retransmitted packets.

Last sections have presented the proposed mechanism and have evaluated its performance in terms of the number of retransmissions. This result has a direct consequence on the bandwidth occupancy.

The bandwidth management and the transmission delays are probably two of the main problems of satellite communications. Indeed, the bandwidth-over-number of potential users ratio is largely lower than for wired Internet or even for wireless LAN. So, the bandwidth must be carefully managed and all possible mechanisms should be used to cope with this problem. A perfect illustration of this is that the next generation of protocols for satellite communications such as e.g. DVB-S2 21 integrate adaptive mechanisms to switch between optimal classes of error correcting codes/modulation.

Another consequence of the reduction of the number of retransmissions is the reduction of the average transmission delay. Clearly, avoiding one retransmission decreases the delay by at least one Round Trip Time (RTT). This RTT being roughly 500 ms for geostationary satellites, any reduction of the number of retransmissions could be very useful for satellite communications. Note that even if large RTTs make it difficult to implement retransmission strategies on satellite communications, several works have shown that ARQ or hybrid ARQ mechanisms can be successfully used in this context [START_REF] Lin | A Hybrid ARQ Scheme with Parity Retransmission for Error Control of Satellite Channels[END_REF] . [START_REF] Azizan | Link Level Performance Analysis of Satellite High Speed Downlink Packet Access (S-HSDPA)[END_REF] Moreover, the intrinsic adaptability of HARQ mechanisms has lead WLAN physical layers designers to integrate them in most of recent standards. [START_REF] Gpp | Technical Specification Group Radio Access Network; High Speed Downlink Packet Access (HSDPA)[END_REF][START_REF]IEEE Standard for Local and metropolitan area networks Part 16: Air Interface for Fixed and Mobile Broadband Wireless Access Systems Amendment 2: Physical and Medium Access Control Layers for Combined Fixed and Mobile Operation in Licensed Bands and Corrigendum 1[END_REF] In our context, we consider an architecture with several users connected to the Internet through a geostationary satellite. Several bidirectional point-to-point communications simultaneously share the same broadcast channel. We consider that the traffic is asymmetric (most of the trafic is sent from the Internet toward the users) with a high rate forward link between the gateway and the user terminals and a low rate return link between the user terminal and the gateway. Note that the return link can also be terrestrial.

In this context, the proposed mechanism can be implemented as described in Section III between the gateway and the user terminals, the gateway acting as the sender.

This proposition applies to many broadcast contexts and some adaptations must be done in order to integrate the satellite-only context. The first point is to determine the layer in which the mechanism can be implemented. At the physical layer, it seems better to consider an error channel where the redundant data are generated with an error correcting code. It can be noted that in the classical mode, a terminal decodes all the sent data (including the data for the other terminals) and then filters them in order to only keep its own data. The modification needed for our mechanism is simply to store temporarily the decoded data to use them for decoding. Thus, it does not imply additional decoding operations upon reception. The implementation of an hybrid-ARQ mechanism was studied [START_REF] Azizan | Link Level Performance Analysis of Satellite High Speed Downlink Packet Access (S-HSDPA)[END_REF] in order to use HSDPA in satellite communications. From this framework, very few modifications are needed to implement our proposal in this context.

To study the implementation of our proposal in DVB-S2, [START_REF] Etsi | Digital Video Broadcasting (DVB); Second generation framing structure, channel coding and modulation systems for Broadcasting, Interactive[END_REF] we can observe that the FECFRAMEs, i.e. the codewords, are generated from information data belonging to different flows. Thus, the implementation of our mechanism can be done simply by implementing classical HARQ type II in DVB-S2. This can be simply done by adding a unique identifier to each FECFRAME header and by allowing some FECFRAME to only contain redundancy which is built from information data sent in a previous FECFRAME.

When implemented at the link layer, the mechanism will be preferably used in erasure mode. The redundant data will then be generated with an erasure code. In a practical context where the link layer is MPE, ULE or GSE, an unique identifier must be added to the data unit headers. Note that this is also a requirement for the implementation of a secure version of these protocols. [START_REF] Noisternig | Lightweight security extension for the Unidirectional Lightweight Encapsulation (ULE) protocol[END_REF] Under this assumption, the generation of repair data unit should be easily implemented since it was already standardized with MPE-FEC 25 or at least proposed. [START_REF] Lei | MPE/ULE-FEC vs GSE-FEC Efficiency Comparison of IP Datagram Transmission over DVB-S2[END_REF]

B. Satellite and terrestrial repeaters

Present or future architectures for TV or radio broadcast from satellite to mobile receivers include terrestrial repeaters (or "gap fillers") that ensure a good reception of the signal in urban areas. Examples of such architectures are MBSAT in South Korea and Japan, Sirius and XM systems in North America or DVB-SH [START_REF] Etsi | Framing Structure, channel coding and modulation for Satellite Services to Handheld devices (SH) below 3 GHz[END_REF] in Europe.

In most of these systems, both satellite and terrestrial links are unidirectional, i.e. there is no return link from the mobile receivers to the satellite or the repeater.

However, with the convergence of technologies (e.g. with 3G technologies), bi-directional links between mobiles and repeaters will probably be available in next-generation architectures.

Under this assumption, we propose the following application of our mechanism: the sender role is split between the satellite and the repeaters.

• the satellite broadcasts the data (toward the mobiles and the repeaters)

• the mobiles send their list of lost data to the repeater • the repeater aggregates the lists, then generates and sends redundant data

• the mobiles recover their lost data and acknowledge the final recovery

The main interest of this scheme is to improve the role of the repeaters in the system. Indeed, in a classical scheme, the repeaters simply broadcasts the data without taken into account the properties of the each mobile. In the propose mode, they can concentrate its efforts to help specific set of mobiles. In addition, the delay between the repeater and a terminal if much smaller than with the satellite, so the overall delay is greatly decreased, thus reducing the impact of retransmissions on upper layers such as TCP.

If we consider the recent standard DVB-SH, [START_REF] Etsi | Framing Structure, channel coding and modulation for Satellite Services to Handheld devices (SH) below 3 GHz[END_REF] it seems easier to implement our mechanism at the link layer instead of the physical layer. Indeed, at the link layer, the standard describes MPE-IFEC, a framework allowing to implement erasure codes. Thus, with the addition of an unique identifier to the MPE data unit, our mechanism should be implemented without any strong modification.

C. Satellite and retransmissions done by peers

In more future architectures, packets lost on the satellite link may be recovered by an ad-hoc network composed of the satellite receivers. In this context, each receiver broadcasts the list of its lost packets and the receivers with good signal reception generate and transmit redundant packets built from these lists of lost packets. This ad-hoc recovery system is quite similar to the network coding COPE system introduced in. [START_REF] Katti | XORs in the air: practical wireless network coding[END_REF] It has the advantage of not requiring any repeater and not having the satellite or gateway manage retransmissions. However such a system should be carefully designed to take into account isolated terminals as well as ad-hoc routing delays.

VII. Conclusion

In this paper we have presented a hybrid ARQ type II mechanism that opportunistically builds a specific Reed-Solomon code based on a minimal set of information packets for each R-receiver transmission block, and produces repair packets using that code. An analytical evaluation of the mean number of retransmissions has been carried out as well as numerical evaluations, that showed that the mechanisms performs very well compared to classical ARQ retransmissions. The mechanism may allow up to a 95% decrease of the amount of retransmitted packets in the case of an error channel. This in turn increases the overall goodput, while the end-to-end delay and jitter are decreased.

This mechanism may be used on any erasure or error broadcast channel, and we have illustrated how to use it in a variety of satellite communications contexts. The most promising one is a mixed satellite/repeater context such as DVB-SH, where the management of the Reed-Solomon code and the transmission or repair packets is done by the repeater, thus even more diminishing the problem of additional end-to-end delays due to retransmissions.

Future work in that direction may include the use of weaker assumptions in our analysis, in particular the introduction of loss correlations in time (e.g. a Gilbert-Eliott error model) or among groups of receivers.

Figure 1 .

 1 Figure 1. Our modified Hybrid ARQ type II mechanism.

Figure 2 .

 2 Figure 2. Retransmissions with erasure and error correcting codes

Figure 3 .

 3 Figure 3. Impact of N and R on the gain

Figure 4 .

 4 Figure 4. Error Channel under-performance wrt Erasure Channel