Sanda Lefteriu
email: sanda.lefteriu@imt-lille-douai.fr

Stefano Grivet-Talocia
email: stefano.grivet@polito.it

Topological Fitting: broadband modeling of passive components via augmented equivalent circuit models

This paper presents a gray-box modeling technique for obtaining physical models of electrical components. These models are obtained through topological fitting, thus yielding equivalent electrical circuit models which are guaranteed to be stable and passive. The algorithm finds the values of the circuit elements given a description of the circuit such that the resulting model fits the frequency response measurements provided for the component. Moreover, in case the circuit is incomplete, the algorithm extends it by automatically generating all series and parallel augmentations and re-optimizing the values of the elements. 68 pt 0.944 in 24 mm 40 pt 0.556 in 14.1 mm 40 pt 0.556 in 14.1 mm 68 pt 0.944 in 24 mm

INTRODUCTION

The increase in the operating frequency of integrated circuits require broadband component models to be easily included in a simulator. Hence equivalent circuit models (ECMs) are designed to take into account parasitic, skin and proximity effects, which are present in high-frequency.

Broadband models are typically developed for each component using physical insight: spiral inductors (Cao et al. (2003); [START_REF] Wang | A wide-band 2 π equivalent-circuit model for spiral inductors on silicon[END_REF]), ceramic capacitors [START_REF] Sullivan | Improved distributed model for capacitors in high-performance packages[END_REF]), thin-film resistors [START_REF] Li | A high-frequency scalable model for thin-film resistor on gan substrate[END_REF]), etc. For large circuits, the number of possible topologies increases dramatically. The choice of element types and values is another issue. Alternatively, frequency responses (Z-,Y-or S-parameters) are used to fit a rational model via blackbox system identification techniques, followed by circuit synthesis [START_REF] Zhaoqing | SPICE Compatible Macromodels of Linear Multiport Networks[END_REF]), however synthesis using only passive elements is not trivial. Passivity is crucial in circuit simulation because a stable but non-passive network may become unstable for certain terminations. This paper aims to make a step towards a fully-automated circuit generation for components described by frequency domain measurements. Oftentimes, one already has some preliminary idea of the topology of an incomplete ECM of a given component or subsystem (nominal circuit). This model can be enhanced by the addition of elements to improve its accuracy. We propose to generate all possible augmentations of the nominal circuit by inserting elements in series or in parallel to the existing ones. To find the optimal value of the added elements, we solve the resulting nonlinear optimization problem via a Sanathanan-Koerner iteration combined with an alternating least-squares.

Previous approaches for automated broadband model generation were not able to ensure stability and passivity by construction. In [START_REF] Kolstad | A new circuit augmentation method for modeling of interconnects and passive components[END_REF], the error between the available measurements and the response of the nominal circuit is modeled by a small number of pole-zero elements, so the final model consists of the ECM together with the augmentation network described as pole-zero elements, which are found by black-box fitting. Stability of the augmentation can be enforced by pole-flipping, while passivity is enforced a-posteriori. [START_REF] Paul | Broadband modeling of high-frequency microwave devices[END_REF] solves the same nonlinear optimization problem as in Eq. ( 15) by expressing the unknown as a rational function and the coefficients of numerator and denominator polynomials are found via least squares. Stability can be enforced by constraining the coefficients to stay positive, while passivity is enforced via perturbation of the purely imaginary eigenvalues of the Hamiltonian matrix pencil [START_REF] Grivet-Talocia | Passive Macromodeling: Theory and Applications[END_REF]). In our approach, we optimize the value of each added circuit element independently in an alternating least-squares setting, thus guaranteeing stability and passivity of the final model.

The paper is structured as follows. Section 2 presents a short review of the modified nodal analysis used to represent electrical circuits. In Sect. 3, we propose to perform all possible series and parallel augmentations of the nominal circuit. Sect. 4 discusses the proposed approach for topological fitting. Lastly, Sect. 5 presents numerical examples validating our approach and Sect. 6 concludes the paper.

CIRCUIT MODELING VIA MODIFIED NODAL

ANALYSIS (MNA) -SHORT REVIEW MNA [START_REF] Ho | The modified nodal approach to network analysis[END_REF]; [START_REF] Vlach | Computer methods for circuit analysis and design[END_REF]) is a powerful technique for automatically formulating (and solving) circuit equations. Multi-port circuits are typically given in terms of their impedance (Z-), admittance (Y-), or scattering (S-) formulation:

• for the impedance formulation, the circuit is currentcontrolled at all ports and the corresponding voltages are measured as outputs:

V = ZI, where V = [ v 1 . . . v p ] T and I = [ i 1 . . . i p ]
T are the vectors of port voltages and currents, respectively, with p, the Margin requirements for the other pages Paper size this page A4 number of ports (we assume, for each port, one input and one output variable).

• for the admittance formulation, the circuit is voltage-controlled at all ports and the corresponding currents are measured as outputs: I = YV. • for the scattering formulation, at each port, the input is the incident wave a = v+R0i 2 √ R0 and the output is the reflected wave b = v-R0i 2 √ R0 , where R 0 is the port reference impedance (typically 50Ω).

Nodes are defined such that there exists a single element between two nodes. For the simple RLC circuit in Fig. 1, nodes have been indicated by circled numbers. Node 0 is the ground (reference for node voltages).

L 1 C 2 R 3 1 2 0 Fig. 1. RLC circuit
The topology of the circuit is described by the incidence matrix: the number of rows is the number of nodes in the circuit and each column corresponds to one branch which contains a single element. Each column is mostly zero, except for the rows given by the two nodes between which this element is connected, which contain 1 and -1.

For the RLC circuit, the incidence matrix is:

A F =   1 2 3 1 1 0 0 2 -1 1 1 0 0 -1 -1   . (1)
Kirchhoff's current law at each node yields the product between the incidence matrix and the branch currents being 0: A F I b = 0. The last equation is the negative of the sum of the previous ones, so we can disregard it and delete the last row to obtain the reduced incidence matrix:

A = 1 2 3 1 1 0 0 2 -1 1 1 . (2)
The Kirchhoff's voltage law expresses branch voltages V b in terms of nodal voltages E: V b = A T E. KCL and KVL describe topological equations, which, together with the constitutive equations, yield the MNA formulation.

We distinguish between the following types of branches:

(1) shunt dynamic: a capacitor (branch 2), (2) shunt resistive: a conductance (no such branch is present in our example), (3) series dynamic: an inductance (branch 1) (4) series resistive: a resistor (branch 3).

The columns of the incidence matrix are partitioned according to the type of element:

shunt dynamic i C = [2], where i C stores their indices, shunt resistive i G = [ ], series dynamic i L = [1] and series resistive i R = [3]: A C = 0 1 , A G = [ ] , A L = 1 -1 , A R = 0 1 . (3)
We also define diagonal matrices storing the values of the elements depending on the type:

C b = [C 2 ] , G b = [ ] , L b = [L 1 ] , R b = [R 3 ] . (4)
For voltage controlled (shunt) elements, the constitutive equations are i G = G b v G and i C = sCv C , with s, the Laplace variable, while for current controlled (series) elements, they are

v R = R b i R and v L = sL b i L .
The unknowns of the MNA system are the node voltages and the branch currents for the series elements. These are included as state variables. In our case:

x(t) = v 1 v 2 i L1 .
The RLC circuit in Fig. 1 has one single port defined between nodes 1 and 0 , corresponding to one input and one output variable. We distinguish between the three formulations:

(1) impedance case: inputs are port currents and outputs are port voltages. Correspondingly, the input and output mapping matrices are

C = [1 0 0 0] , B = C T .
(2) admittance case: inputs are port voltages and outputs are port currents. Since the input current is not present in the vector of generalized state variables, we add a series resistive element (with zero-valued resistance), whose current will be our input variable. Correspondingly, the matrices after the update are:

A R = 0 -1 1 0 , R b = R 3 0 .
The input and output mapping matrices are

C = [0 0 0 0 1] , B = C T .
(3) scattering case: inputs are the incident waves and outputs are the reflected waves. Incident waves can be represented by series Thevenin sources with internal resistance R 0 , hence the incidence matrix A R should be updated with such new elements, and the diagonal R should be updated with a block containing the reference R 0 on the diagonal:

A R = 0 -1 1 0 , R b = R 3 R 0 .
The input and output mapping matrices are

B = 0 0 0 0 2 √ R 0 , C = 1 2 √ R0 0 0 0 - √ R0 2
.

The descriptor matrices of the MNA formulation are

E dae =   A C C b A T C 0 0 0 L b 0 0 0 0   , A dae =-   A G G b A T G A L A R -A T L 0 0 -A T R 0 R b   (5)
and the transfer function is given by the usual formula

H(s) = C (sE dae -A dae ) -1 B, (6)
which can represent the Z-,Y-or S-parameter formulation. H(s) is a rational matrix function of variable s and size p. This dynamical system is described by a differential algebraic equation with an invertible matrix pencil (A dae , E dae ), a singular E dae and the output equation having D = 0: Margin requirements for the other pages Paper size this page A4

E dae ẋ(t) = A dae x(t) + Bu(t), (7)
y(t) = Cx(t). (8

CIRCUIT AUGMENTATION

Given a nominal (starting) circuit, it is possible to insert elements in series or in parallel according to these rules:

(1) parallel augmentations with shunt elements between any two existing nodes if no such element is already connected to these nodes, (2) series augmentations of existing branches with series elements if no such element is already present.

Parallel augmentation of a given circuit described by its incidence matrix and types/values of branch elements is simply performed by adding an element connected to existing nodes. The new element is described by a new column appended in the incidence matrix, its type and value.

On the other hand, series augmentation is performed by adding an element in series to an existing branch. For this, the branch is split into two series-connected branches by adding a node appended to the list of existing nodes. The original branch is then reconnected between one of the original nodes and the new node, and the new branch is connected between the new node and the second original node. Consequently, the row size of the incidence matrix grows by one and some entries are modified to account for the reconnection. An additional column is appended to the incidence matrix to account for the new element.

To minimize user interaction, we propose to generate all possible series and parallel augmentations of a nominal circuit according to the two rules above. This operation can be performed once or several times, however, the complexity increases significantly after each step. The elements added by augmentation will account for parasitics, skin or proximity effects in high-frequency and will yield an automatically-generated broadband ECM.

ELEMENT PERTURBATION

We are provided with frequency responses of components obtained by direct measurement or using an electromagnetic field solver. Using minimal user insight, a nominal circuit is defined. This circuit is then augmented following Sect. 3. Lastly, we perturb the value of each circuit element to match the response of the new system to the given measurements.

Single element

We treat the case of a single element first and we distinguish between two scenarios based on its type:

• shunt element: (5) shows that the stamp of the element in the MNA system is recovered through the appropriate column of the incidence matrix. If the element is connected between nodes k and l, the selector vector ξ indicates its connectivity

ξ =         0 1 ← k th entry 0 -1 ← l th entry. 0 (9)
If node k or l is the ground node, the corresponding row is deleted, hence ξ contains only 1 or -1.

• series element: (5) shows that the value of the element directly multiplies the corresponding current unknown in the state-space vector, hence

ξ = 0 n e L 0 nr or ξ = 0 n 0 n l e R (10)
for perturbing an inductance or a resistance, respectively. 0 n is a null column vector of size n, the number of nodes, 0 nr , of size n r , the number of series resistive elements and 0 n l , of size n l , the number of series dynamic elements. Moreover, e L is a unit column vector with 1 at the position of the corresponding inductor and similarly for e R .

Denoting sE dae -A dae as G, the transfer function is H(s) = CG -1 B (note that G contains the frequency dependency). The perturbation induced in the element value induces a perturbation in G as follows:

G = G 0 + ∆αξξ T , (11)
where G 0 is the matrix containing the nominal values and ∆α is the perturbation in the value of the element. If we treat a dynamic element, a capacitance or an inductance, then ∆α = s∆C i or ∆α = s∆L i . For static elements, ∆α = ∆G i or ∆α = ∆R i , respectively.

The Sherman-Morrison formula [START_REF] Golub | Matrix computations[END_REF]) allows to express G -1 :

G -1 = G -1 0 - G -1 0 ∆αξξ T G -1 0 1 + ∆αξ T G -1 0 ξ . (12)
Hence, the transfer function is

H(s) = CG -1 B = CG -1 0 B H0(s) -∆α Fc(s) CG -1 0 ξ F b (s) ξ T G -1 0 B 1 + ∆α ξ T G -1 0 ξ K(s) = H 0 (s) -∆α F(s) 1 + ∆α • K(s) , ( 13
) with F(s) = F c (s)F b (s) (note that K(s) is scalar).
We wish to minimize the deviation between the measurements provided, denoted by H i , and the transfer function of the perturbed system evaluated at the measurement frequencies ω i , where ω i = 2πf i , for i = 1, . . . , N , with N , the number of samples provided:

E = N i=1 H(iω i ) -H i 2
F , so, after substitution of ( 13), the error becomes

E = N i=1 H 0 (iω i ) -∆α F(iω i ) 1 + ∆α • K(iω i ) -H i 2 F = N i=1 ∆H(iω i ) -∆α F(iω i ) 1 + ∆α • K(iω i ) 2 F , with ∆H i =H 0 (iω i )-H i . We redefine F(iω i ) as F i ∈ C p×p and K(iω i ) as K i ∈ C, so the error expression becomes E = N i=1 ∆H i -F i ∆α 1 + ∆α • K i 2 F . ( 14
)
The optimization problem is simply

min ∆α E = min ∆α N i=1 ∆H i -F i ∆α 1 + ∆α • K i 2 F . ( 15
)
The optimization parameter ∆α appears in the numerator and denominator, making this problem nonlinear. Employing the Sanathanan-Koerner (SK) iteration (Sanathanan Margin requirements for the other pages Paper size this page A4

and Koerner ( 1963)), we relax the problem by using ∆α from the previous iteration in the denominator, thus yielding a linear problem which can be solved to determine the current ∆α:

min ∆α N i=1 ∆H i + (∆H i K i -F i ) ∆α k 1 + ∆α k-1 • K i 2 F (16) via least-squares: A LS ∆α k ≈ b LS , with A LS =       vec ( (X 1 )) vec ( (X 1 )) . . . vec ( (X N )) vec ( (X N ))       , b LS =       vec ( (Y 1 )) vec ( (Y 1 )) . . . vec ( (Y N )) vec ( (Y N ))       , ( 17
)
where vec is the vectorization of the p × p matrix in the argument, , are the real and imaginary parts and lastly,

X i = ∆HiKi-Fi 1+∆α k-1 Ki and Y i = -∆Hi 1+∆α k-1 Ki .

Several elements

For perturbing several elements, we employ the alternating least squares technique. Each element is optimized independently by keeping the other parameters fixed and applying Sect. 4.1. This iterative process is run until convergence or a maximum number of iterations is reached.

The pseudocode is detailed in Algorithm 1.

NUMERICAL EXAMPLES

RLC circuit with parasitic elements

The first example is a theoretical example (with unrealistic values) in which the goal is to identify the parasitic elements of an inductor and a capacitor in an RLC circuit.

The circuit is shown in Fig. 2, with the parasitics indicated in red, the parasitic elements having values much smaller than those of the circuit elements. We model the inductor with value 1.5H as having an equivalent series resistance of .05Ω. Similarly, the capacitor with value .22F is modeled together with a series resistance of .005Ω and a series inductance of .03H.

1.5H .05Ω .22F .005Ω .03H 60Ω - + V in i out

Fig. 2. RLC circuit with parasitics

We consider admittance parameters for the circuit in Fig. 2, by having the circuit voltage-controlled and measuring the current in the first branch. We generate 200 measurements of the admittance parameters for this circuit between ω = 1 and ω = 10rad/s. We start from the RLC circuit with no parasitic elements, which we augment using all possible series augmentations by R and L elements and parallel augmentations by G and C elements. Afterwards, topological fitting is applied to the augmented circuit by perturbing the values of the augmented elements and the values found are indicated next to each element in Fig. 3 Algorithm 1: top fit Input : (f i , H i ), ∀i = 1, . . . , N : measurements, R 0 : diagonal matrix with input impedance at each port: ∞ (impedance), 0 (admittance), real positive (scattering), A: reduced incidence matrix, nodes: structure with node information, list el: list of elements as a cell array of strings {'R','L','C','G'}, val el: vector of nominal values, ind el: vector of indices for the elements to be perturbed opts: structure with tolerance and maximum iterations number Output: new val el: vector of new numerical values E dae , A dae : updated descriptor-form

1 [A C , A G , A L , A R , C b , G b , L b , R b , i C , i G , i R , i L ] =
get top(A,list el,val el) // Obtain circuit topology as in (3) and (4)

2 nx = #rows(A)+size(L b )+size(R b ) 3 [A R , R b , B, C] = add top ports(nx,nodes,R 0 ,A R ,R b )
// Based on Z-/Y-/S-formulation, extend A R and R as needed and create B and

C 4 [E dae , A dae ] = construct DAE(A C ,A G ,A L ,A R ,C b ,L b , 5 G b ,R b ) // Constructs A and E as in (5) 6 ξ ξ ξ = create xi(ind el,list el,i C ,i G ,i L ,i R ,A C ,A G ,A L ,A R ,
size(E dae ))// selector vectors as in (9),(10) 7 nIter 1 = 0, 1 = ∞// outer loop 8 α α α = 0 // initial guess 9 ∆α α α = 0// stores increments 10 α α α =val el(ind el)// stores current values 11 while nIter 1 <opts.maxIter & 1 >opts.tol do (the conductance is represented as a rectangle, while the resistor is displayed using the zig-zag symbol). Lastly, Fig. 4 shows the Y-parameter measurements of the circuit in Fig. 2, which, together with the RLC circuit with no parasitics, are used as input parameters to our approach to recover the circuit in Fig. 3.

= L norm L, C = C norm Ĉ, R = R norm R.
We first fit a π configuration (Fig. 6) to the first half of the frequency band to find the elements matching this lowfrequency behaviour. The initial guesses for the elements' values are 10 -7 (this is the nominal circuit for this step). The values found with top fit yield a model which matches accurately the first half of the measurements (shown in Fig. 7). In the next step, we duplicate this circuit to obtain a configuration, thus accounting for proximity effects (Cao et al. (2003); [START_REF] Wang | A wide-band 2 π equivalent-circuit model for spiral inductors on silicon[END_REF]). We distribute the inductor value found in the first step and use the value found in Fig. 6 divided by two as initial guess for the two inductors, while the initial guess for the last capacitor is 10 -7 . This constitutes the nominal circuit for this step. Topological fitting is applied again on the whole frequency band. The circuit found is depicted in Fig. 8. The measurements, the response of the nominal model and the response corresponding to the circuit in Fig. 8 are shown in Fig. 9. This circuit fits the measurements quite accurately, with a slight discrepancy still observed in the imaginary part of S 12 .

  This work was supported by the ACRI foundation through the Young Investigator Training Program award and the TD COST Action TD1307 European Model Reduction Network (EU-MORNET) as a Short Term Scientific Mission.

  Fig. 3. Augmented RLC circuit with the values found

Fig. 4 .Fig. 5 .

 45 Fig. 4. Results for the RLC circuit with parasitics Comparable results were obtained when adding white Gaussian noise with 60dB signal-to-noise ratio to the measurements. The circuit found is shown in Fig.5while results are similar to those in Fig.4. The noise level was chosen such that the parasitic elements with values in the 10 -3 range are not embedded in noise.

Fig. 7 .

 7 Fig. 6. π-configuration found from the low-frequency behavior of the spiral inductor (normalized component values)

  Fig. 8. 2π-configuration for the spiral inductor (normalized component values) In the last step, we generate all possible series and parallel augmentations of the 2π circuit of Fig. 8 (this is the

  )

	68 pt
	0.944 in
	24 mm
	40 pt	40 pt
	0.556 in	0.556 in
	14.1 mm	14.1 mm
	68 pt
	0.944 in
	24 mm

9.2 • 10 -6

1.8 • 10 -5 10 -5

Fig. 10. Circuit obtained after applying topological fitting on the circuit obtained by all possible series and parallel augmentations of the 2π configuration (normalized component values)

CONCLUSION

This paper aims to make a step towards automated generation of equivalent circuit models for components. These models are valid in a broad-band and are stable and passive by construction. The algorithm works well for simple topologies with few circuit elements, but may run into difficulties if the size grows too large or if the topology of the nominal circuit is not suitable for the application at hand. In the future, we consider combining the SK process with a Newton iteration and use regularization for alternate least-squares to improve and speed up convergence.

Cao, Y., Groves, R.A., Huang, X., Zamdmer, N.D., Plouchart, J.O., Wachnik, R.A., King, T.J., and Hu, C.


Laurie André

Stéphane Abanades
email: stephane.abanades@promes.cnrs.fr

Laurent Cassayre

Experimental and Thermodynamic Study of Co-Fe and Mn-Fe Based Mixed Metal Oxides for Thermochemical Energy Storage Application

Metal oxides are potential materials for thermochemical heat storage, and among them, cobalt oxide and manganese oxide are attracting attention. Furthermore, studies on mixed oxides are ongoing, as the synthesis of mixed oxides could be a way to answer the drawbacks of pure metal oxides, such as slow reaction kinetics, loss-in-capacity over cycles or sintering, selected for thermochemical heat storage application. The addition of iron oxide is under investigation and the obtained results are presented. This work proposes a comparison of thermodynamic modelling with experimental data in order to identify the impact of iron oxide addition to cobalt oxide and manganese oxide. Fe addition decreased the redox activity and energy storage capacity of Co 3 O 4 , whereas the cycling stability of Mn 2 O 3 was significantly improved with added Fe amounts above 20 mol% while the energy storage capacity was unchanged. The thermodynamic modelling method to predict the behavior of the Mn-Fe-O and Co-Fe-O systems was validated, and the possibility to identify other mixed oxides becomes conceivable, by enabling the selection of transition metals additives for metal oxides destined for thermochemical energy storage applications.

INTRODUCTION

The purpose of thermal energy storage in solar power plants is to store concentrated solar heat in the form of chemical energy during on-sun hours and then to use it during off-sun hours, thus enabling energy production according to needs. Thermal energy storage (TES) can be achieved using three main routes: latent, sensible and thermochemical heat storage [START_REF] Abedin | The open renewable energy journal[END_REF]. This study focuses on thermochemical heat storage, which shows advantages over latent and sensible heat storage, including higher energy storage densities, possible heat storage at room temperature in the form of stable solid materials, and long term storage in a large temperature range (400-1200°C) with a constant restitution temperature. Concentrated solar power (CSP) can provide the heat source required for endothermal/exothermal reversible reactions involved in thermochemical energy storage (Eq. 1, 2) in order to store solar energy in the form of chemical bonds. When combined with CSP systems the continuous production of electricity in power plants becomes possible. The reaction enthalpy is stored in the reaction products during the heat charge (Eq. 1), and can be released by reversing the reaction during the discharge (Eq. 2).

MO n + Heat

 MO n-1 + ½ O 2 (g) (1)
MO n-1 + ½ O 2 (g)  MO n + H e a t

Cobalt oxide, Co 3 O 4 , and manganese oxide, Mn 2 O 3 , are two of the most studied metal oxides considered as promising materials for thermochemical energy storage. Experiments showed that Co 3 O 4 is the most suited raw material given the fast reaction kinetics and complete reaction reversibility [START_REF] Agrafiotis | [END_REF][3][4][5][6]. However, the cost and potential toxicity of cobalt oxide requires the development of other alternative materials. Optimization of materials reactivity is required for the other metal oxide species by using e.g. doping strategies, controlled synthesis techniques for tailored morphology or stabilization with inert materials to alleviate sintering effects. The effect of the morphology of the material has been studied, for instance, by Agrafiotis et al. (2015) [START_REF] Agrafiotis | [END_REF]3] who achieved thirty redox cycles with porous foam of Co 3 O 4 . Karagiannakis et al. (2014) [4] studied cobalt oxide, Co 3 O 4 /CoO, in the form of flowthrough pellets and stated that these structures show better reaction kinetics than cobalt oxide as a powder, due to better heat transfer characteristics of the structured material. They also studied manganese oxide, Mn 2 O 3 /Mn 3 O 4, in the form of flow-through pellets. The material showed better reaction kinetics than those obtained with its powder form. Because of the very slow re-oxidation kinetics of manganese oxide, the re-oxidation occurs generally in two steps, one during the cooling and the next one during the re-heating of the next cycle [4]. The addition of a secondary metal oxide can be a way to modify the properties of a pure oxide.

Improvement of materials properties and flexibility such as reaction kinetics and reaction temperature tuning can be obtained by the addition of dopants. For example, Mn 2 O 3 is known for its slow re-oxidation kinetics [4], which affects its cycling stability. The doping with transition metals can relieve this issue [3,7]. The present study aims at investigating the effect of Fe addition on the performances and thermodynamic properties of Co and Mn-based mixed metal oxides which have been identified as interesting candidates for thermochemical energy storage. The optimal mixed metal oxide composition to reach the highest performances for TES is addressed, encompassing high reaction enthalpy, stability upon cycling, low temperature gap between reduction and re-oxidation, and reaction in a temperature range adapted to CSP. The addition of transition metals is investigated as a way to answer the drawbacks of pure metal oxides, such as slow reaction kinetics [4], sintering [8] or loss-in-capacity over cycles [9]. Block et al. (2014) [7] reported that both addition of iron oxide to cobalt oxide or addition of cobalt oxide to iron oxide results in lower enthalpies of reaction compared to those of the pure oxides. However, they estimated that cobalt oxide doped with 10% of iron oxide possesses higher reduction/oxidation reversibility than pure cobalt oxide and still shows high enthalpy of reaction. Pagkoura et al. (2014) [10] showed that cobalt oxide with 10-20w% of iron oxide shows good thermo-mechanical stability over ten redox cycles. Carillo et al. (2014) [8] showed that the addition of iron oxide does not allow avoiding the sintering encountered with Mn 2 O 3 but helps to increase the heat storage density of the material. Furthermore, the addition of iron helped to stabilize and enhance the oxidation rate of manganese oxide over an experiment of thirty redox cycles. According to their study, the fastest and most stable oxidation reaction was obtained for Mn 2 O 3 doped with 20% Fe.

The possibility to predict the optimal amount of transition metal to be used and its impact on the material's properties, such as the transition temperatures, the phases formed and encountered during the heating, the structure changes and the oxygen exchange capacity, is required to assess the suitability of the redox process. For a better understanding and prediction of the behavior of Co-Fe-O and Mn-Fe-O systems at high temperatures, this work proposes a comparison of thermodynamic modelling performed with FactSage and of experimental studies for validating thermodynamic predictions. An advantage of using metal oxides as thermochemical energy storage material is the use of air as the heat transfer fluid that can be processed in an open-loop system. For this reason, the thermodynamic modelling calculations take into account an atmosphere of PO 2 = 0.21 atm. The experiments are conducted in the same atmosphere, using 20%O 2 /Ar as a representative gaseous mixture. The mass loss (O 2 storage capacity), the temperatures of reduction and oxidation and the enthalpies of reaction were measured and compared to calculations given by thermodynamics.

SYNTHESIS AND CHARACTERIZATION

All the powders were synthesized via a modified Pechini method, using metal nitrates, citric acid and ethylene glycol in aqueous solution [7]. The Fe molar ratio presented in percentage in the following is based on a Fe/(Fe+Co) and Fe/(Fe+Mn) mole ratio. The powders were all calcined at 750°C for 4h in air in order to fully eliminate the organic fraction and residues of the synthesis, and to stabilize the structure. The complete removal of organics was confirmed by the absence of extra mass loss during the first heating step of thermogravimetric analyses. Higher calcination temperatures were not considered to avoid the reduction of the materials during the synthesis. The powders were characterized by X-ray diffraction (XRD) before being studied in redox cycles. XRD analysis was performed (Fig. 1) at room temperature using a PANalytical XPert Pro diffractometer (CuKα radiation, λ = 0.15418 nm). X-ray diffraction measurements of θ-θ symmetrical scans were made over an angular range of 10 to 80°. The step size and the time per step were fixed at 0.01° and 20 s respectively. The contribution from AlKα 2 was removed and the X-ray diffractograms were recorded and studied using the PANalytical software. The instrumental function was determined using a reference material (SRM 660, lanthanum hexaboride, LaB6 polycrystalline sample) and can be expressed by a polynomial function [11]. Pure Co 3 O 4 and Mn 2 O 3 , and Co 3 O 4 and Mn 2 O 3 with iron incorporated into the structure, were synthesized. The more the iron added to Co 3 O 4 , the more amorphous the material gets (Fig. 1.b). For both Co 3 O 4 and Mn 2 O 3 , the peaks for iron oxide are covered by noise or by the base-oxide peaks, and are only detected for the material containing the largest amount of iron. The synthesized oxide phases were corresponding to the structures of (i) a cubic spinel solution (Co 1-x Fe x ) 3 O 4 (Fig. 1a andb) and (ii) a bixbyite solution (Mn 1-x Fe x ) 2 O 3 [START_REF] Lai | Mixed iron-manganese oxide nanoparticles[END_REF] (Fig. 1c andd). Traces of iron oxide (Fe 3 O 4 and Fe 2 O 3 , respectively) were detected for the samples with the highest Fe content (Fig. 1b andd).

EXPERIMENTAL RESULTS

Experimental data concerning transition temperatures during redox process, oxygen storage capacity of redox materials and reaction enthalpies were obtained by using simultaneous thermogravimetric analysis (TGA) and differential scanning calorimetry (DSC), using a Netzsch STA 449 F3 System. Reduction-oxidation cycles were performed in a 20%O 2 /Ar atmosphere, between 800°C and 1050°C, with pure cobalt oxide and cobalt oxide mixed with 5 mol%, 10 mol% and 25 mol% Fe (Fig. 3), and between 750°C and 1050°C with pure Mn 2 O 3 and manganese oxide mixed with 10, 20, 30, 40 and 50 mol% Fe. Each TGA profile also contains a separate reduction step under Ar, in order to get information about the reduction temperature of the material under inert atmosphere. The heating rate was 20 °C/min for the reduction step, and the cooling rate was 10°C/min for the re-oxidation step. The gas flow rates used were 10 Nml/min for O 2 and 40 Nml/min for Ar. About 50 mg of powder was placed into an alumina crucible and then analyzed in TGA coupled with DSC. A series of three cycles was performed for each studied formulation as it proves sufficient in order to determine their suitability for achieving reversible reactions and to assess a possible chemical deactivation of the materials as a result of a loss of redox stability during thermal treatment.

Experimental data were compared with equilibrium calculations in order to validate the thermodynamic modelling approach. The FactSage software was combined with the FToxid database, which includes a full description of the thermochemical properties of the oxide phases in the Fe-Mn-O and Fe-Co-O systems, in order to calculate the behavior of mixed oxides intended for thermochemical energy storage applications. The theoretical mass loss, Δm, was calculated at atmospheric pressure and under air between 800°C and 1050°C, and the enthalpy of phase transition was determined.

Cobalt Oxide

Cobalt oxide is considered as a promising thermochemical energy storage material due to its high gravimetric storage density and high cycling stability when compared with other metal oxide candidates. The reduction temperature of Co 3 O 4 to CoO in air atmosphere was measured between 880 and 930°C [4,6,10,13]. A section of the Co-Fe-O phase diagram, presented in Fig. 2a, has been calculated for PO 2 = 0.21 atm. The phase diagram shows that the iron content greatly influences the composition and the amount of the monoxide phase formed when heating the spinel phase. As an example, Fig. 2b presents the influence of the temperature on the conversion rate for cobalt oxide mixed with 20 mol% Fe: a temperature of about 1380°C is required to fully convert the spinel phase into the monoxide phase. 4 with addition of various amounts of iron (Fig. 3). Co 3 O 4 containing the highest amount of iron (Fig. 3a) shows the lowest reduction extent (the mass loss, Δm, during heating up to 1050°C is the lowest). Increasing the amount of iron in the material thus reduces the mass of O 2 it is able to release. As observed in the phase diagram (Fig. 2a), the material composition transitions from a spinel phase to a spinel plus monoxide phase when heating above ~900-1000°C for Fe content in the range considered in this study (5-25%). The higher the iron amount added to the material, the lower the reduction to cobalt monoxide and the lower the amount of O 2 released by the material. In thermochemical energy storage, the enthalpy of the reaction is related to the oxygen storage capacity of the material. When the redox activity of the material is lowered, there is a decrease of reaction enthalpy. Thus, with more iron added, the gravimetric energy storage density is decreased compared with pure Co 3 O 4 , resulting in less storage capacity per gram of material. The theoretical estimations of Δm obtained from thermodynamics show a linear decrease with increasing Fe content and are perfectly in agreement with the obtained experimental values (Fig. 4a). However, the cycling stability of cobalt oxide suffers from the addition of iron oxide, as the re-oxidation conversion rate is slightly decreasing over multiple cycles for Co-Fe mixed oxides, which is not the case for pure The influence of the incorporation of iron oxide on reaction temperatures was also studied and temperatures at peak reaction rate, i.e. the temperatures at which the reaction kinetics reached its maximum, were measured. A gradual increase of the reduction and oxidation temperatures is observed with increasing amount of iron in the material (Fig. 4b). An effect of Fe addition on the temperature gap between reduction and oxidation is also visible. As both temperatures rise with more Fe addition, it can be noticed that, at 25 mol% of added Fe, the gap between the reduction and oxidation temperatures has been reduced (Fig. 4b). The reduction of this gap in temperature is interesting for large-scale application since it reduces the amount of energy spent for the heating and cooling of the system. The measured enthalpies of reaction are decreasing with increasing amount of Fe added, which agrees with the decrease in Δm. The highest enthalpy measured is 471.6 kJ/kg for pure Co 3 O 4 , and the lowest is 194.4 kJ/kg for Co 3 O 4 with 25 mol% Fe, which correlates well with the observed Δm differences in the TGA measurements of the two formulations (Fig. 4a); while, for pure Co 3 O 4 , the theoretical enthalpy is about 844 kJ/kg according to thermodynamics [14], and about 747 kJ/kg according to FactSage models. The noteworthy discrepancy between the measured and theoretical enthalpies is the result of Co

Manganese Oxide

The Mn 2 O 3 /Mn 3 O 4 couple was studied by TGA coupled with DSC with and without addition of iron. The reduction step of Mn 2 O 3 was observed in the range of 920-1000°C and the notably slow re-oxidation was observed in the range of 850-500°C with a gravimetric energy storage density of 110 kJ/kg. The re-oxidation happens in two steps, the first one being during the cooling and in between 700°C-500°C and the second one being during the reheating and in the range of 500-850°C [10]. The temperature of the TGA program used here is ranging between 750°C and 1050°C. The calculated Mn-Fe-O phase diagram is presented in Fig. 5, for PO 2 = 0.21 atm. It can be observed that, when heated under air, manganese oxides mixed with less than 20% of iron go through two consecutive two-phase zones before reaching full reduction into a single spinel phase. In contrast, above about 20 mol% Fe, the material fully transforms from the bixbyite Mn 2 O 3 phase to the cubic spinel Mn 3 O 4. Thus the amount of O 2 exchanged, Δm, remains the same whatever the amount of iron, until reaching 50 mol% (limit studied here) (Fig. 5a). This is well illustrated with the example of the phase transition of Mn 2 O 3 with 20 mol% Fe in Fig. 5b. The phase transition is fast during heating and takes place promptly in a very narrow range of temperature, at about 960°C, according to experimental results. TGA shows that the sample with the composition going through the two-phase zones, (i.e. containing 10 mol% Fe added to Mn 2 O 3 ), is unable to regain its full mass during oxidation, as the Δm lost during the first reduction step cannot be recovered during the re-oxidation step. Inversely, the other compositions, from 20 mol% Fe to 50 mol% Fe, are capable of regaining their lost O 2 mass in a complete reversible way. They all share the same Δm, as predicted by the Mn-Fe-O phase diagram. In addition, the measured enthalpy of reaction is in agreement with theoretical estimations, for pure Mn This tendency is consistent with the experimental results, showing that the maximal values of Δm are similar to the theoretical ones. Mn 2 O 3 shows cycling stability issues due to sintering. This phenomenon is reflected by the average experimental Δm, showing a decrease of redox activity during cycling. Mn 2 O 3 mixed with 10 mol% Fe also loses cycling stability, regaining only 31% of its lost mass. However, from 20 mol%Fe to 50 mol% Fe added, a great improvement is observed, since the average Δm value remains close to the maximum Δm value and to the theoretical value (Fig. 6a.). The addition of iron to Mn 2 O 3 reduces the sintering of the material and enhances the cycling stability, as already observed previously [8].

As for the influence of Fe incorporation in Mn 2 O 3 on the reaction temperatures, a similar tendency as with cobalt oxide is obtained. The reduction and oxidation temperatures of Mn 2 O 3 increase with increasing amount of added iron (Fig. 6b). The addition of iron oxide to Mn 2 O 3 enhances its cycling stability (in the range 20-50 mol% Fe) and also changes its reaction temperature. This way, the reaction temperature can be tuned. In addition to an increase of the temperature, further addition of iron also reduces the gap in temperature between the reduction and the oxidation step.

In summary, this study showed that the addition of iron to manganese oxide reduces the sintering of Mn 2 O 3 and enhances the cycling stability of the material. Also, the higher the amount of iron added to Mn 2 O 3 , the higher the reaction temperatures. In addition, the gap in temperature between the reduction and the oxidation decreases with higher iron content.

CONCLUSION

The effect of Fe addition in Co and Mn-based oxides was studied experimentally and results were compared with thermodynamic calculations. The addition of iron to Mn 2 O 3 is promising since it can be used to tune the temperature of redox reactions, to reduce the gap in temperature hysteresis between the reduction and the oxidation step, as well as to enhance the re-oxidation kinetics and cycling stability of the material by countering the sintering issue of Mn 2 O 3 . Conversely, the incorporation of iron to Co 3 O 4 has negative effect on the redox performances since the maximum amount of O 2 exchanged is lowered when the amount of iron added is increased. The increasing amount of iron added to Co 3 O 4 reduces the maximum oxygen mass exchanged during a redox cycle. The addition of iron to Co and Mn-based oxides also results in an increase of the reaction temperatures, while slightly lowering the gap in temperature between the reduction and oxidation step. The comparison of the obtained experimental data with FactSage thermodynamic modelling proves that this method is effective for predicting the behavior of mixed oxides since it results in consistent phase transition temperatures and reliable oxygen storage capacities. It could be further used for selecting other transition metals to enhance the properties of mixed metal oxides for thermochemical energy storage application.

FIGURE 1 .

 1 FIGURE 1. XRD analysis of (a) Co 3 O 4 with 5 mol% Fe, (b) Co 3 O 4 with 25 mol% Fe, (c) Mn 2 O 3 with 10 mol% Fe and (d) Mn 2 O 3 with 50 mol% Fe.

FIGURE 2 .

 2 FIGURE 2. (a) Co-Fe-O phase diagram under air, (b) example of phase composition and conversion rate evolution with the temperature for Co 3 O 4 with addition of 20 mol% Fe. The TGA of Co 3 O 4 shows the variations in the amount of O 2 exchanged, Δm, for Co 3 O 4 with addition of various amounts of iron (Fig. 3). Co 3 O 4 containing the highest amount of iron (Fig. 3a) shows the lowest reduction extent

  Co 3 O 4 (Fig. 4a.). In Fig. 4a are represented, along with calculated values, the maximum values of Δm measured for each composition, which means the highest amount of O 2 released during reduction, and the average Δm for the three redox cycles. For pure Co 3 O 4 , the average Δm remains the same as the theoretical value and the maximum value measured experimentally.This means that the O 2 released and regained did not change, confirming the material cycling stability. For the samples containing iron, since the average Δm falls below the theoretical and maximum value, a slight loss in O 2 exchange capacity upon cycling is evidenced (Fig.4a).

FIGURE 3 .

 3 FIGURE 3. TG analysis of Co 3 O 4 /CoO + (a) 25 mol% Fe, (b) 10 mol% Fe, (c) 5 mol% Fe

  3 O 4 structural changes attributed to spin-state change [14]. A similar decrease in enthalpy with increasing Fe content was observed by Block et al (2016) [7] as they measured an enthalpy of 433 kJ/kg for Co 3 O 4 with 13 mol% Fe, and a lower enthalpy of 165 kJ/kg for Co 3 O 4 with 33 mol% Fe.

FIGURE 4 . 4 /

 44 FIGURE 4. (a) Average and maximal experimental Δm (%) compared to theoretical Δm between 800°C and 1050°C (under atmosphere condition of 20% O 2 ), (b) Experimental temperatures at peak reaction rate for pure cobalt oxide and mixed Co-Fe oxides.

FIGURE 5 .

 5 FIGURE 5. (a) Mn-Fe-O phase diagram under air, (b) example of phase composition evolution with the temperature for Mn 2 O 3

2 O 3 ,

 23 of 190.1 kJ/kg[15], and about 207.78 kJ/kg according to FactSage calculations. An average of 185.2 kJ/kg is measured for the samples with compositions between 20 mol% Fe and 50 mol% Fe. The samples of pure Mn 2 O 3 and with 10 mol% Fe show similar enthalpy, 145 kJ/kg and 144.2 kJ/kg respectively, measured for the first reduction step with full conversion. The highest enthalpy of 193.3 kJ/kg is measured for Mn 2 O 3 with 20 mol% Fe. Experimental and theoretical mass variations, Δm, were compared. It is observed that Δm remains stable with the addition of iron to manganese oxide, in agreement with the phase diagram.

FIGURE 6 .

 6 FIGURE 6. (a) Average and maximal experimental Δm (%) compared to theoretical Δm (under atmosphere condition of 20% O 2 ), (b) Experimental temperatures at peak reaction rate for pure manganese oxide and mixed Mn-Fe oxides.


Mauricio Iturralde
email: mauricio.iturralde@imt-atlantique.fr

Thomas Galezowski
email: thomas.galezowski@societedugrandparis.fr

Xavier Lagrange
email: xavier.lagrange@imt-atlantique.fr

Performance of Mobile Relays in Loaded Conditions for Railway Transportation

Keywords: QoS, LTE, 4G, mobile relay

Deploying mobile relays in public transportation is a simple yet effective way to avoid the electromagnetic insulation within vehicles and to increase the Quality of Service (QoS) perceived by passengers, which is sometimes low. Mobile relaying can be done by using an LTE/EPC 100%-compatible architecture. However, this solution induces extra-overhead. In this paper we evaluate the QoS of an LTE mobile relay architecture for public railway transport systems for two representative services: client-server requests and voice communications. We compare the performance of a direct transmission against a mobile relay architecture for different types of requests and different load conditions. This work, therefore, evaluates the mobile relay performance in terms of load time, throughput, packet loss ratio and end-to-end latency. Our findings show that a mobile relay architecture highly improves the QoS performance for train passengers. Furthermore, the gain is greater as the load increases.

I. INTRODUCTION

The use of wireless broadband services is rising significantly with the deployment of Long Term Evolution (LTE) networks and the generalisation of smart phones, tablet computers and other new mobile devices. People make intensive use of these devices when they are on public transport vehicles such as buses, trams, or trains. It is anticipated that, by the year 2020, it will be fairly common to have up to 50 active vehicular User Equipments (UEs) per bus and up to 300 active vehicular UEs per train [START_REF] Popovski | Scenarios, requirements and kpis for 5g mobile and wireless systems. The METIS project: Mobile and wireless communications Enablers for the Twenty-twenty Information Society[END_REF].

The service quality in public transport is far from satisfactory. Vehicles are usually well shielded with coated windows, which leads to a rather high penetration loss between outdoor and in-vehicle. The penetration loss can be as high as 25 dB, and even goes up to 35 dB in the Shanghai high-speed magnetic levitation train [START_REF] Liu | Position-based modeling for wireless channel on highspeed railway under a viaduct at 2.35 ghz[END_REF]. Traditionally, the UEs inside the public transport system are connected to macro base stations via wireless links, in which the penetration loss severely attenuates the signal quality and decreases the achievable data rate. Deploying mobile relays with both an outdoor antenna to communicate with the cellular network and an indoor antenna to provide a good coverage to passengers is a natural solution to provide a high QoS.

Société du Grand Paris is in charge of designing and constructing the new 200km of fully automated metro lines around Paris, connecting the 3 Paris airports as well as the main suburban and innovation areas. These "Grand Paris Express" lines, interconnected with the existing Paris subway, regional trains and bus lines, will carry around 2 Million passengers per day. Société du Grand Paris is seeking to provide continuous high quality telecommunication services to passengers in Grand Paris Express stations and inside trains. Thus, Société du Grand Paris is interested in developing new technologies such as mobile relays in order to serve this purpose.

In [START_REF] Kerdoncuff | Mobile relay for lte: proof of concept and performance measurements[END_REF] we showed with a testbed using real radio transmissions that a mobile relay architecture can be easily implemented with standard Evolved Packet Core (EPC) and with full compatibility with 3GPP recommendations. However, since this solution induces extra overhead and extra latency, the QoS of this solution needs to be evaluated.

The objective of this paper is to evaluate the performance of the mobile relay architecture where a train is loaded with several users using HyperText Transfer Protocol (HTTP) and Voice over Internet Protocol (VoIP) services.

The remainder of this paper is organised as follows. In Section II we present the state-of-the-art of mobile relay in LTE. In Section III we describe the architecture and the protocol stack when mobile relays are used. In Section IV we analyse the advantages and limits of mobile relays. In Section V we describe our simulation scenario within the traffic models and we also present our QoS parameters. In Section VI we present and analyse our findings and finally Section VII concludes our document.

II. RELATED WORKS

Previous works have focussed on the architectural aspect of Mobile Relay [START_REF] Krendzel | Lte-a mobile relay handling: Architecture aspects[END_REF] [5] [START_REF] Chen | Analysis and improvement of mobility procedures for mobile relays in lte networks[END_REF], data-rate gain and QoS aspects. In [START_REF] Krendzel | Lte-a mobile relay handling: Architecture aspects[END_REF] the authors outline and sum up the main network architecture aspects for mobile relay handling. They analyse two types of mobile relay architecture. This paper takes into account the handover aspects as well. In [START_REF] Chen | Mobile relay in lte-advanced systems[END_REF], the architecture for mobile relay is presented. The key techniques of supporting mobile relay are investigated, such as the group mobility, the local service support, the multi Radio Access Technology (RAT) and Radio Access Network (RAN) sharing, with the corresponding solutions. In [START_REF] Chen | Analysis and improvement of mobility procedures for mobile relays in lte networks[END_REF] authors propose to keep the same protocol stack as for fixed relays but introduce the concept of a global tunnel, which gathers several tunnels for the mobile relay architecture. This architecture is compatible with the GPRS Tunneling Protocol (GTP). The authors show that without adding any protocol in the stack, it is possible to group several handovers while keeping the possibility of managing the tunnel of each UE individually.

Other works have pointed their attention to the data rate gain aspect [START_REF] Chen | Downlink data rate gain provided by a mobile relay for lte-advanced[END_REF], [START_REF] Sui | The potential of moving relays-a performance analysis[END_REF], [START_REF] Chen | Performance analysis of cellular networks with mobile relays under different modes[END_REF], [START_REF] Tang | Coverage performance of joint transmission for moving relay enabled cellular networks in dense urban scenarios[END_REF] and [START_REF] Chen | Downlink capacity gain analysis of mobile relay in lte-advanced network[END_REF]. The authors analyse in [START_REF] Chen | Downlink data rate gain provided by a mobile relay for lte-advanced[END_REF] the data rate gain provided by a mobile relay deployed in public transportation systems. Also an analysis of the penetration loss factor between outdoor and in-vehicle is presented. They found that the loss factor has an important influence on the data rate gain when using a relay architecture. In [START_REF] Chen | Downlink capacity gain analysis of mobile relay in lte-advanced network[END_REF] the authors analyse the data rate gain when using mobile relays in public transportation. The loss factor and the on-board user's throughput gain when avoiding it is also studied. They found that avoiding the train loss factor, the on-board users receive a 30% higher bit rate. In [START_REF] Sui | The potential of moving relays-a performance analysis[END_REF], authors studied the performance of dual-hop transmission with a mobile relay and with a fixed relay by means of a theoretical analysis and simulations and proved the benefit of mobile relays when the penetration loss is above 25 dB. In [START_REF] Chen | Performance analysis of cellular networks with mobile relays under different modes[END_REF], stochastic analysis is used to compute the capacity of a cell when mobile relays are used. In [START_REF] Tang | Coverage performance of joint transmission for moving relay enabled cellular networks in dense urban scenarios[END_REF] joint transmission by both a mobile relay and a macro-cell base station is proposed to improve the data rate of non-vehicular UEs and the outage probability is derived.

The QoS aspect is studied in [START_REF] Calle-Sanchez | Theoretical analysis and performance simulation of in-band lte mobile relays in railway environments[END_REF]. Authors studied the mobile relay capacity improvement for on-board train users and its impact on the overall network performance. This paper analyses the gain in terms of throughput when transmitting video streaming for only 10 on-board train users. Some other studies consider mobility protocol issues. In [START_REF] Chen | Mobile relay in lte-advanced systems[END_REF], several architectures based on the 3GPP analysis [START_REF]3rd Generation Partnership Project. Evolved Universal Terrestrial Radio Access (E-UTRA); Study on mobile relay[END_REF] are described. In [START_REF] Meng-Shiuan Pan | An enhanced handover scheme for mobile relays in lte-a high-speed rail networks[END_REF], [START_REF] Huang | Mobile relay based fast handover scheme in highspeed mobile environment[END_REF], the authors optimized the handover for mobile relays deployed on high speed trains by using the predictability of movement on a railway. In [START_REF] Pietsch | Moving Relays and Mobility aspects[END_REF], the authors discussed several architecture alternatives for mobile relays and proposed intra mobility area mobile relay handover and inter mobility area mobile relay handover procedures. The handover latency of these procedures are analysed. However, this proposal requires extra planning effort on the operator side and some upgrades of the standard 3GPP based protocols.

The Third Generation Partnership Project (3GPP) investigated several possible architectures for mobile relays in [START_REF]3rd Generation Partnership Project. Evolved Universal Terrestrial Radio Access (E-UTRA); Study on mobile relay[END_REF]. Some solutions require the modification of several protocols and thus an important standardisation effort. However, the first alternative (called alt1) described in [START_REF]3rd Generation Partnership Project. Evolved Universal Terrestrial Radio Access (E-UTRA); Study on mobile relay[END_REF] is based on the fact that nodes of the EPC such as Mobile eNB (eNBm), Mobile Management Entity (MME) just need Internet Protocol (IP) connectivity and that IP connectivity can be simply provided by an LTE/EPC network. In [START_REF] Kerdoncuff | Mobile relay for lte: proof of concept and performance measurements[END_REF] we proved this concept with real radio transmissions.

To our knowledge, there has been as yet no systematic examination of interactive and VoIP services in mobile relays. The question remains, however, how those services perform under a mobile relay architecture in loaded conditions. The wireless access for passengers is provided by the use of two nested LTE/EPC networks as shown in Figure 1. The first one is called the Track Evolved Packet Core (EPCt) and acts as a backhaul network. Nodes of this network are identified by suffix t (eNBt, MMEt, SGWt, etc.). We name a terminal of this network a Transport Operator terminal Equipment (TOE). There is typically one TOE in every train. The second network is the general network. It manages the UEs of passengers and could be owned by operators. It has eNBm as well as standard eNB. Note that eNBms have access to the general network through the track network: an eNBm is connected to a TOE and can thus exchange messages with Global Mobile Management Entity (MMEg) and Global Serving Gateway (SGWg).

It is worth mentioning that all these elements are standard LTE components, and that the use of this architecture requires no adaptation of protocol stacks or hardware.

There can be several general networks that use the same LTE/EPC track network in order to allow all mobile operators of a country to provide the service on-board. There can be several eNBms (one per operator) but Radio Access Network (RAN) sharing can also be used to avoid a duplication of embedded hardware in trains. Note that Figure 1 shows a functional architecture. The operator of the track network can be a mobile operator and several entities can be physically in the same node (e.g. MMEt and MMEg are in the same piece of equipment)

The protocol stack when a mobile relay is used is shown in Figure 2 on the radio interface and on the S1-interface. Due to the use of two LTE/EPC networks there is an additional level of encapsulation and thus extra-headers compared to a standard architecture. The following headers are added to each packet: 8-12 bytes for GTP [18], 8 bytes for the User Datagram Protocol (UDP) and 20 bytes for the Internet Protocol (IP) (with IPv4). The GTP header length is 12 only when packet numbering is activated in GTP. We consider the worst case and thus assume a 12-byte header, which gives a 40-byte total overhead. As the length of the internal packets is variable, an important question is to determine whether this 40-byte overhead is negligible or not compared to the user data average size for popular services such as telephony and HTTP. In this section, we analyse the possible benefits brought by the deployment mobile relays. We consider two representative services, namely client-server requests and telephony.

A. Interactive services

We consider a generic interactive service as shown in figure 3 with a request-response protocol. The UE sends a request to a server. The server then transfers a file, which can be a web page, a video or a text. The application protocol is HTTP and the transport protocol is the Transmission Control Protocol (TCP).

Fig. 3: HTTP client-server model

The main QoS indicator is the loading time defined as the total time to transfer the whole file from the time the request was sent. To get a better insight into the procedure, we break down the loading time into the server response time and the file reception time (see figure 3).

Due to the additional transmission on the backhaul compared to a direct configuration, there is one more radio hop with a mobile relay. The server response time is thus increased. According to [START_REF] Pries | An http web traffic model based on the top one million visited web pages[END_REF], a website's size can be from 150 KB up to 8 MB. Thus IP packets are set to their maximum length, which is generally 1500 bytes. The 40-byte extra overhead is thus negligible. The file reception time can be dramatically reduced thanks to a much better radio link budget.

B. Telephony service

In 4G networks, telephony is based on VoIP. End users are generally sensitive to the packet loss rate and the one-way delay. ITU-T G.114 [20] recommends a maximum of 150 ms one-way latency in order to avoid any problems due to possible echoes and to keep excellent interactivity between the two parties. Since this includes the entire voice path, part of which may be on the public Internet, the LTE/EPC network should have a transit latency less than 150 ms. As for the interactive service, using a relay increases the one-way delay.

The relative cost of the extra-overhead on the backhaul can be not negligible. Though there are several possible modes for the speech codec, the most representative one is the GSM Enhanced Full Rate (EFR) codec: A 244-bit data block is sent every 20 ms, which gives a rate of 12.2 kbps. There is generally a cyclic redundancy check (CRC) that is added. The payload is thus 32 bytes. The sizes of the header added by the VoIP protocol stack are 12 bytes, 8 bytes and 20 bytes for RTP, UDP and IPv4, respectively. We assume that the PDPC, RLC and MAC layers add headers whose total length is 6 bytes [START_REF]ns-3 simulator[END_REF]. In standard VoIP over LTE, the required bitrate on the radio interface is (6 + 40 + 32) × 8/20 = 31.2 kbps. With the relay extra overhead (see figure 4), the required bitrate on the backhaul is (6 + 40 + 40 + 32) × 8/20 = 47.2 kbps. The increase is 47.2/31.2, that is +51%. The question is to know whether the throughput increase really compensates the additional overhead.

A. Simulation Environment

For performing our tests we used network simulator version 3 (ns-3) [START_REF]ns-3 simulator[END_REF]. This is a discrete-event network simulator for Internet systems that has a well developed LTE Module for The LTE Model includes the LTE Radio Protocol stack (RRC, PDCP, RLC, MAC, PHY) and the EPC Model includes core network interfaces, protocols and entities. These entities and protocols reside within the SGW, PGW and MME nodes, and partially within the eNB nodes. We use the standard NS3 LTE model with direct transmission between UEs in the train and the network as the reference mode.

The LTE module has been modified to implement the mobile relay architecture scenario (see Figure 1). In particular, the encapsulation and decapsulation mechanisms are fully simulated. We tested the performance of the mobile relay architecture and we compared them against a direct transmission scenario.

In our scenario, a train moves in a linear way with a velocity of 100 kmph as shown in Figure 5. Base stations (eNBt) are regularly deployed along the track and have a 600-m radius. The intersite distance is thus 1200 m. Handover is perfect both for the direct and the relay modes (no packet loss and no call drop). The eNbt and the eNbm transmit on different frequencies. Therefore, the inter-cell interference is avoided.

Two train loads are considered : 20 and 100 mobile users. We assume 60% of the passengers make HTTP requests and the remaining 40% have VoIP communications. The simulation parameters are presented in detail in Table I.

B. HTTP Traffic Model

We considered different mean sizes for the file : 20 000, 100 000 and 500 000 bytes. Inspired by [START_REF] Pries | An http web traffic model based on the top one million visited web pages[END_REF], we assume the size follows a Weilbull r.v., whose CDF is

F X (x) = 1 -e -(x/λ) k . ( 1
)
where k > 0 is the shape parameter and λ > 0 is the scale parameter of the distribution. The value of k used is 0.8. The value of λ is fixed in order to get the chosen mean size for files. For mean sizes of 20 000, 100 000 and 500 000 bytes the values of λ are 17 652, 88 261 and 441 305, respectively. After each file transfer, the user reads the content. The reading time is modelled by an exponential r.v. with a mean value of 39.7 seconds. At the end of reading time, a new request is sent to the server by the terminal.

C. VoIP Traffic Model

We simulate VoIP communications with an Enhanced Full Rate (EFR) codec or equivalently an AMR-WB codec in 12.65-kpbs mode (see [START_REF]3rd Generation Partnership Project. Universal mobile telecommunications system (umts); lte; codec for enhanced voice services (evs)[END_REF]). The source traffic generates a 32-bit data block every 20-ms. There is no voice activity detection. Hence, the flow is constant bit rate. The RTP/UDP/IP headers are added to each data block. We assume the calls are set up prior to the simulation and did not consider any call arrival and call release. No signalling is simulated.

D. QoS Parameters

For the HTTP flows we define the following parameters:

Server Response Time: This indicator is computed for each request u as

η u = t fp,u -t rq,u (2)
where, t fp,u is the time stamp of the first received packet by the UE and t rq,u is the time stamp of the request sent by the UE to the server (see figure 3). Page Reception Time: This indicator is computed as

ξ u = t lp,u -t fp,u (3)
where t lp,u is the time stamp of the last packet received by the UE. Page Load Time: The page load time is simply the sum of the server response time and the page reception time:

τ u = η u + ξ u (4)
Throughput: This parameter measures the throughput obtained when receiving the HTTP file. This parameter is calculated as:

α = 1 U U u=1 B u t lp,u -t fp,u (5)
where U is the total number of requests during the simulation, B u is the amount of bytes received after request u.

And for the VoIP flows we analyse the latency and the Packet Loss Ratio (PLR).

Packet Loss Ratio: This parameter measures the end-toend packet losses. This parameter is calculated as:

ρ = 1 U U u=1 P tx,u -P rx,u P tx,u (6)
where P tx,u is the number of transmitted packet by the server, P rx,u is the number of received packets at the UE.

VI. RESULTS AND ANALYSIS

In this section we first analyse the QoS performance of HTTP client-server file transfer by regarding several parameters such as: file load time, server response time, and throughput. Secondly we look at the VoIP performance based on QoS parameters such as: PLR and end-to-end latency.

Our figures represent the statistical distribution of values by using box plots with the upper quartile, median and lower quartile. The whiskers represent the maximum and minimum values excluding outliers. Outliers are the points which fall more than 1.5 times the interquartile range above the third quartile or below the first quartile.

A. Results of HTTP flows

Figure 6 depicts the effect of the mobile relay in terms of page load time when the train is loaded by 60 users downloading HTTP files (within a total of 100 users). For 500-kB files, the median load time is reduced from 7 s to 2.5 s and is thus divided by 3. The maximum load time is reduced by 30 s. Users are generally very sensitive to such a reduction of the load time. Even for moderate-size file for which the load time is by nature small, it is reduced by the use of a relay. In all cases, the maximum load time is always smaller in the relay mode (1 s vs 2 s for 20-kB files, 4 s vs 15 s for 100-kB files, 22.5 vs 52 s for 500-kB files).

The average sever response time is increased by about 20 ms when a relay is used as shown in figure 7. This is due to the supplementary radio hop. However, the response time can be larger than 80-ms (and even sometimes larger than 120 ms) in the direct mode due to bad radio conditions as opposed to the relay-mode where the maximum value is about 75 ms. This means that using a relay is beneficial even for clientserver applications in which a very small amount of data is downloaded.

Figure 8 shows the throughput when downloading an HTTP file. When the train is close to the base station, radio conditions are excellent and the rate is high. Thus, the maximum throughput in the direct mode is close to the one in the relay 3 Mbps for 100-kB files). However, the minimum throughput is lower in direct mode at about 0.1 Mbps in comparison to relay mode at 0.25 Mbps. Even without specific QoS policy, the mobile relay helps to offer a minimum bit rate (of course except in highly saturated conditions with a huge number of passengers). The average throughput is much larger with the relay mode. The throughput gain increases as the file size is larger. Due to the slow-start mechanism, the TCP protocol takes greater advantage of the best channel conditions reached by the relay architecture when the file is larger.

Figure 9 shows the effect of the mobile relay in terms of page load time when the train is loaded by 12, 60 and 120 users downloading HTTP files of 100-kB. For 12 users the median load time is less than 1 s for both the direct and the relay modes. For such low load, the relay does not provide any benefit but does not degrade the QoS.

The maximum load time is reduced by the use of a relay from 16 s to 4 s for 60 users and from 39 s to 13 s for 120 The median is divided by 2 and the maximum load time is divided by 3 for 120 users. For moderate to high loads, there is a clear interest to use mobile relays. Note that a train loaded by 120 users in relay mode shows a shorter load time than a train loaded by 60 users in direct mode. This means that the use of a mobile relay allows serving about twice the number of users with a slightly higher QoS.

B. Results of VoIP flows

Figure 10 shows the distribution of the latency. Since in the relay architecture packets have more hops, it increases the latency from about 25 ms to 25-37 ms. As the latency is less than 40 ms, this guarantees a good voice calling interaction since the latency threshold is 150 ms. Furthermore, a simple Proportional-Fairness scheduler is used. It does not give any priority to VoIP packets. By using a dedicated bearer for VoIP packets on the backhaul, it is possible to give a higher priority to voice packets and to limit the latency increase. The PLR for VoIP flows was measured for each simulation. It is not shown in this paper because only very few packets were lost and using boxplots for displaying results has no statistical significance. For both direct and relay modes, the loss is less than 0.01%. This confirms that the radio conditions of the direct mode are not so bad and our comparison for all services is fair: the configuration parameters we chose represents a reference case (i.e. the direct mode) with a good quality coverage scenario. In such conditions, the gain for interactive services is clearly visible. It can be larger if the coverage inside the train is poor in the direct mode.

VII. CONCLUSIONS

This study has focused attention on the analysis of the QoS performance of mobile relays in a standard LTE architecture for a railway scenario. We evaluated the QoS of HTTP and VoIP services when a train is loaded up to two hundred users simultaneously using the LTE Network.

Implementing a mobile relay solution in trains highly improves the quality of interactive services and keeps an excellent quality for VoIP services.

It is clear from the current study that avoiding penetration losses improves the QoS for the passengers, even if the mobile relay architecture generates an additional overhead of 40 bytes to each transmitted packet. The gain for interactive services is more important for long downloaded files and for high load conditions.

Since the QoS in a relay architecture depends of the behavior of both the mobile eNb and the eNBs covering the track, future works should focus on scheduling techniques for real-time and non-real time services. Also the performance of handovers should be studied.

Fig. 1 :

 1 Fig. 1: Architecture of a network with a mobile relay

Fig. 2 :

 2 Fig. 2: Protocol header in a mobile relay

Fig. 4 :

 4 Fig. 4: VoIP packet

Fig. 5 :

 5 Fig. 5: Mobility Scenario

Fig. 6 :Fig. 7 :

 67 Fig. 6: Boxplot of http load time for different file sizes

Fig. 8 :

 8 Fig. 8: Boxplot of throughput for http flows

Fig. 9 :

 9 Fig. 9: Boxplot of http load time for different loads

Fig. 10 :

 10 Fig. 10: Boxplot of VoIP latency

TABLE I :

 I Simulation Parameters


Qipeng Song

Xavier Lagrange

Loutfi Nuaymi

Performance of selection combining macro diversity with outage probability constraint in LPWAN

Keywords: macro diversity, selective combining, LPWAN

In low power wide area networks (LPWAN), a packet transmitted by a device can be received and decoded by all surrounding base stations (BS). If at least one of the base stations decodes the packet, the latter is delivered to the system. This scheme is called selection combining (SC) macro diversity and has been studied in the literature. However, the existing research only cares about the average packet-loss probability over the entire networks. In this paper, we consider as the main performance indicator the outage probability, which better takes into account the real availability of the service in any part of the area: the outage probability. Given an outage probability constrain, we study the minimum required BS spatial density of SC macro diversity based LPWAN. The numerical results show that if both target packet loss rate and outage probability are 10%, the minimum BS spatial density required by SC macro diversity is at most half of that required by traditional BS attach method. The performance gain is more significant when the network load increases.

Introduction

Low Power Wide Area Network (LPWAN) is regarded as a promising solution to handle future machine type communication (MTC) traffic [START_REF] Song | Survey of radio resource management issues and proposals for energy-efficient cellular networks that will cover billions of machines[END_REF]. There are several technologies, but most representative ones (Sigfox and LoRaWAN) are based on Aloha because of its simplicity [START_REF] Raza | Low power wide area networks: An overview[END_REF].

In cellular networks, a device always attaches to the Base Station(BS) for which the received power averaged over all fading realizations is the strongest and then transmits packets to its attached BS. This scheme is referred to as Best BS attach mode. However, a device in LPWAN works in a different way: it sends a packet in broadcast mode without BS attach procedure and benefits from macro reception diversity, which is defined as the capacity of several BS to receive the same packet (see Fig. 1). Each BS autonomously and independently decodes the packets and then sends the decoded packets to the core network. That is to say, each BS is equipped with a packet decode function, and the backhaul link between BS and core network transmits only when a packet is received at the BS. The core network is in charge of duplicate received packets removal (e.g., by comparing the identity and message content conveyed in packets). A packet is successfully delivered if at least one BS decodes the packets. This scheme is referred to as selection-combination-based macro diversity, simply written as SC macro diversity . It is currently used by Sigfox and LoRaWAN [START_REF] Raza | Low power wide area networks: An overview[END_REF]. In the literature, the performance of ALOHA is usually evaluated with indicator (e.g. packet loss rate) that is averaged over the entire networks. This kind of indicator gives a global view about the quality of service in LPWAN networks but still have its limitations. For example, with a target network level packet loss rate such as 10%, it is very possible that devices in the neighborhood of BS never fail while devices at the border of BS coverage area have packet loss rate higher than 10%.

In this paper, we consider a performance indicator that better takes into account the real availability of the service in any part of the area: the outage probability . This indicator is defined as the probability that at any location a device has a link-level packet loss rate greater than a target value. For example, the statement that a LPWAN network has an outage probability of 10% with target packet loss rate 1%, means that at any location, the probability that a device has a packet loss rate higher than 1% is 10%. With outage probability, we can study the minimum BS spatial density if a LPWAN operator needs to guarantee a certain level of grade of service (i.e. outage probability). These studies can help LPWAN operators to deploy their network in an efficient way.

Combining techniques were studied to improve the performance of cellular networks in the 1990s. For example, the performance of SC macro diversity is studied for GSM networks with a remote antennas system by using hexagonal topology [START_REF] Jung | Exploitation of intracell macrodiversity in mobile radio systems by deployment of remote antennas[END_REF]. The packet loss rate and system capacity in ALOHA networks are usually analyzed with Stochastic geometry. The performance of SC macro diversity for BS with multiple receiver branches was studied in [START_REF] Pena-Martin | Performance of Selection Combining Diversity in η -µ Fading Channels With Integer Values of µ[END_REF]. Stochastic geometry has served as a powerful tool to study the performance of LPWAN. In this field, Baccelli et al. [START_REF] Baccelli | An Aloha protocol for multihop mobile wireless networks[END_REF] use Laplace transform for interference analysis for slotted ALOHA. Haenggi et al. [START_REF] Haenggi | Interference in large wireless networks[END_REF] extensively study the outage probability in SINR-based capture model for slotted ALOHA. B laszczyszyn et al. [START_REF] B Laszczyszyn | Stochastic analysis of non-slotted ALOHA in wireless ad-hoc networks[END_REF] open the door to study pure ALOHA with stochastic geometry. In our previous work [START_REF] Song | Evaluation of Macro Diversity Gain in Long Range ALOHA Networks[END_REF], we neglected the background noise and obtained closed-form expression of average packet loss rate in SC macro diversity based Aloha LPWAN and showed the performance gain compared to best BS attach. However, to our best knowledge, SC macro reception diversity with outage probability for ALOHA systems has not been studied in literature.

The remaining of this paper is organized as follows: Sec. 2 introduces the system model in which we take into account the Rayleigh fading, shadowing, background noise and capture effect. Sec. 3 analyzes the successful packet transmission probability for a given pair of device and BS based on some stochastic geometry research results. These analyzes cover both pure and slotted ALOHA. Sec. 4 presents our analysis about the minimum required BS spatial density under outage probability constraint. Sec. 5 presents numerical results and discussions. Sec. 6 concludes this paper.

System model

For any random variable (RV) X, let f X (x), F X (x), L X [s], φ X (ω) be its probability density function (PDF), cumulative distribution function (CDF), Laplace transform (LT) and characteristic function (CF), respectively.

Distribution of Nodes and Traffic Model

We consider a large wireless network over a two-dimension infinite plane. The locations of terminals form a stationary Poisson point process (PPP) Φ m = {x i } on the plane R 2 , where x i refers to the coordinates of device with label i. The spatial density of Φ m is λ m . Similarly, the locations of base stations also form a stationary PPP Φ b = {y j } with spatial density λ b .

Slotted ALOHA and Pure ALOHA

In the system model, we consider pure ALOHA (P-ALOHA) and slotted ALOHA (S-ALOHA). For slotted ALOHA, the time domain is equally divided into slots with duration T slot . In each slot, each device independently decides to transmit a packet with probability p. The propagation delay is assumed to be much smaller than T slot . Hence, there is a global slot synchronization over the whole network. In pure ALOHA, devices send packets without synchronization, but we still use parameter p. The packet generation process can be seen as an internal slotted system in each device with p the probability to transmit in a slot of duration T slot . At a given time, the locations of terminals that are transmitting a packet form a thinned PPP with spatial density pλ m . For slotted and pure ALOHA, we define the normalized load (per BS) as L = pλ m /λ b .

Propagation model

The propagation model is based on Okumura-Hata with both shadow and Rayleigh fading. We assume that MTC devices transmit with identical power level denoted by P t . The received power P r of a packet at the BS is given by:

P r = P t Kr g -γ H10 σ dB χ/10 = P t Kr g -γ H exp σ dB ln(10) 10 χ , (1)
where K is a constant which depends on the antenna characteristics and the average channel attenuation, r g refers to the Euclidean distance between the transmitter and the receiver, γ is the path-loss exponent, H is an exponentially distributed RV with mean 1, and χ is a standard normal variable. For the sake of simplification, we define σ = σ dB ln(10)/10. We assume that H and χ are both constant during a packet transmission, and mutually independent for different links.

Capture effect is taken into account to determine whether a BS could decode a received packet. Thus, the transmission success probability p s for a transmitreceiver pair is defined as follows:

p s = P SINR = P r I + σ 2 n ≥ θ T , (2)
where I refers to the suffered cumulative interference during packet transmission, σ 2 n is the background noise power.

Displacement Theorem

To facilitate the analysis, we use the displacement theorem, which is formulated as a lemma in [4, lemma 1].

Lemma 1. For a homogeneous PPP Φ ⊂ R 2 with spatial density λ b , if each point r g ∈ Φ is transformed to r ∈ R 2 such that r = X -1 γ r g , where {X} are i.i.d., such that E X -2 γ < +∞, the new point process Φ ⊂ R 2 defined by the transformed points r is also a homogeneous PPP with density λ b = λ b E X 2 γ .
We thus rewrite (1) as follows:

P r = P t KHr -γ , with r = P t K exp(-χ/γ)r g , (3)
which indicates that the performance analysis in the initial PPP with shadowing is equivalent to that in a transformed PPP of intensity

λ b = λ b E e 2 γ χ = λ b e 2σ 2 γ 2
without shadowing but with the modified distance r. The PDF of r, proved in [START_REF] Andrews | A tractable approach to coverage and rate in cellular networks[END_REF], is given as follows:

f r (r) = 2πλ b exp(-λ b πr 2 )r, r ∈ [0, +∞] , (4)
In the following, all distances are the modified ones in the transformed PPP.

Link-level Transmission Success Probability

We first study the transmission success probability p s (r) for a given uplink between a device with label x 0 and a device with label y 0 . The distance between x 0 and y 0 is denoted by r. For the considered device x 0 , all the other interfering devices constitute a PPP Φ m \ {x 0 }. Combining ( 1) and ( 2), we have:

p s (r) = P P t KHr -γ xj ∈Φm\{x0} P t KH xj r -γ xj + σ 2 n ≥ θ T , (5)
Using Slivnyak's theorem [START_REF] Vaze | Random Wireless Networks: An Information Theoretic Perspective[END_REF][Theorem 2.3.3], (5) can be further simplified:

p s (r) = P P t KHr -γ xj ∈Φm P t KH xj r -γ xj + σ 2 n ≥ θ T . (6)
Let I = xj ∈Φm P t KH xj r -γ xj , which is the cumulative interference suffered by x 0 . As shown in [START_REF] Haenggi | Interference in large wireless networks[END_REF], p s (r) can be expressed in terms of Laplace transform of cumulative interference L I (s) at point θ T r γ :

p s (r) = P H ≥ (I + σ 2 n P t K )θ T r γ = exp(- σ 2 n θ T P t K r γ ) +∞ 0 exp(-yθ T r γ )f I (y)dy = exp(- σ 2 n θ T P t K r γ )L I {θ T r γ } , (7)
Now the problem is down to calculate the Laplace transform of cumulative interference suffered by the BS at the origin. Interference I in 2 is constant during a packet transmission for S-ALOHA, but may vary with respect to time for P-ALOHA because other MTCs start or stop their transmissions during the packet transmission. With advanced transmission techniques (e.g., interleaving, robust channel coding, etc.), I is the average value within the considered packet transmission and the access is called Pure Aloha with average interference (Pa-ALOHA). For both access schemes, the Laplace transmission of I, proved in [START_REF] B Laszczyszyn | Stochastic analysis of non-slotted ALOHA in wireless ad-hoc networks[END_REF], is as follows:

L I (s) = exp -pλ m πA exp 2σ 2 /γ 2 s 2/γ , ( 8
)
where A is defined as follows:

A = Γ(1 -2/γ) Γ(1 + 2/γ) , for S-ALOHA 2γ γ+2 Γ(1 -2/γ) Γ(1 + 2/γ) , for Pa-ALOHA
where Γ (•) is the gamma function.

For P-ALOHA without advanced transmission techniques, I is the maximum interference value during the packet transmission and is called maximum interference (Pm-ALOHA). Its closed-form expression of LT does not exist [START_REF] B Laszczyszyn | Stochastic analysis of non-slotted ALOHA in wireless ad-hoc networks[END_REF]. In [11, Sec.III-B], we propose to approximate I yj by un upper bound that is the sum of interference levels at the start and end time of the considered packet transmission. With our proposed proposed upper bound, ( 8) is extended to cover the case of Pm-ALOHA with

A = 2 Γ(1 -2/γ) Γ(1 + 2/γ).
By combining ( 7) and ( 8), we can generalize the link-level success transmission probability p s (r) for all ALOHA cases (slotted and pure ALOHA) :

p s (r) = exp(- σ 2 n θ T P t K r γ ) exp(-pλ m πAθ 2 γ T e 2σ 2 γ 2 r 2 ) = exp (-η T r γ ) exp -T r 2 (9)
where σ 2 n is the background noise power, θ T is the capture ratio, P t refers to the device transmit power, K is a positive constant related to propagation model and constant A is defined in [START_REF] Lagrange | Réseaux GSM: des principesa la norme[END_REF]. For the simplicity of notation, let η T =

σ 2 n θ T PtK , T = pλ m πAθ 2 γ T e 2σ 2 γ 2 .

Minimum required BS density

The transmission success probability p s (r) paves the way to study the minimum required spatial density λ b,min for best BS attach and SC macro diversity.

Best BS Attach

Consider a device at the origin. If best BS attach method is applied, according to ( 9), the link level packet loss rate p f,b (r) is:

p f,b (r) = 1 -exp (-η T r γ ) exp(-T r 2 ), (10)
where

η T = σ 2 n θ T PtK , T = pλ m πAθ 2 γ T e 2σ 2
γ 2 . Let P f,max be the target network loss rate and P outage,b be the outage probability (the subscript b refers to the best BS attach method). Since the link level packet loss rate p f,b (r) is a function of r, thus, the proportion of devices that have higher packet loss rate than P f,max , namely the outage probability P outage,b , can be expressed as follows:

P outage,b = +∞ 0 1 r [p f,b (r) > P f,max ]f (r)dr = +∞ rc f (r)dr = 1 - rc 0 f (r)dr = exp(-λ b e 2σ 2 γ2 πr c 2 ), (11)
where 1 r [•] is a indicator function about r and its PDF is given in ( 4), r c is the critical distance that can be numerically calculated from equation p f,b = P f,max . The critical distance r c is calculated by substituting p f,b (r c ) = P f,max into [START_REF] Raza | Low power wide area networks: An overview[END_REF].

1 -exp (-η T r γ c ) exp(-T r 2 c ) = P f,max (12)
Although there is no closed-form expression for critical distance r c , the latter can be easily numerically calculated by leveraging the fact that p f,b (r) is a monoincreasing function with respect to r.

With the obtained r c and a given outage probability, the minimum BS spatial density λ b,min can be obtained by inversing [START_REF] Song | Evaluation of Macro Diversity Gain in Long Range ALOHA Networks[END_REF]:

λ b,min = -ln(P outage,b ) 1 π exp(2σ 2 /γ2)r 2 c , (13)
where ln (•) is the natural logarithm operator.

Without background noise

In an interference-limited LPWAN system, the background can be reasonably neglected. Closed-form expression r c can be obtained from ( 12) by letting η T = 0:

r c = ln(1 -P f,max )/ -pλ m πAθ 2 γ T e 2σ 2 γ2
1/2 [START_REF] Vaze | Random Wireless Networks: An Information Theoretic Perspective[END_REF] With substitution of ( 14) into ( 13), λ b,min is as follows:

λ b,min = pλ m ln(P outage,b )Aθ 2/γ T ln(1 -P f,max ) . ( 15
)

SC Macro Reception Diversity

Different from best BS attach method, whose link-level packet loss rate only depends on the distance to the nearest BS, the performance of SC macro reception diversity depends on a distance vector r i , i = 0, 1, 2.... Thus, we need to know the characteristics of such a vector. However, the distance to the nearest BS r 0 has the most impact to SC macro reception diversity. Thus, to reduce the complexity and keep the tractability, for a given device, its packet loss rate can be approximately expressed in terms of r 0 (instead of the distance vector to all surrounding BS r i , i = 1, 2... ).

Let p f,sc (r 0 ) be link-level packet loss rate with respect to r 0 . Recall that in case of SC macro reception diversity,

p f,sc (r 0 ) ≈ (1 -p s (r 0 ))E ri∈Φ b (1 -p s (r i )) , with r i ∈ [r 0 , +∞] , (16)
where r i refers to the distance between the device and BS with label i. Note that r i should be not less than r 0 . The impact of BS other than the nearest one to the packet loss rate is thus reflected by term 16) is actually the Probability Generating FunctionaL (PGFL) of PPP Φ b for BS, which states for some function f (x) that:

E ri∈Φ b (1 -p s (r i )) . Term E ri∈Φ b (1 -p s (r i )) in (
E x∈Φ f (x) = exp(-λ( R 2 (1 -f (x))dx))
Thus, it can be simplified as follows:

E ri∈Φ b (1 -p s (r i )) = exp -πλ b e 2σ 2 γ2 +∞ r0 2 exp(-η T r γ 2 -T r)dr (17)
With substitution of ( 17) into (16), we have:

p f,sc (r 0 ) = (1 -e -η T r γ 0 -T r 2 0 ) exp -πλ b e 2σ 2 γ2 +∞ r0 2 exp(-η T r γ 2 -T r)dr (18)
By letting p f,sc (r 0 ) = P f,max and substituting (13) into (18), the critical distance r c can be numerically calculated due to the fact that p f,sc (r 0 ) is a monoincreasing function (even though its closed-form expression does not exist). Then we can use formula [START_REF] Sornin | LoRaWAN specification[END_REF] to calculate the minimum BS spatial density, given an outage probability and target network packet loss rate.

Numerical result

The system settings to get the numerical results are resumed in Tab. 1. Some explanations are given as follows: to evaluate the impact of the background noise to the system performance, we need to determine the values of P t , K, σ 2 n . For one LPWAN device running on ISM band (e.g. with 868 MHz), from a regulatory point of view, the effective radiated power (ERP) may not exceed 14 dBm (or 25 mW) in any direction [START_REF] Sornin | LoRaWAN specification[END_REF], we thus assume that the transmit power P t is set as 14 dBm.

About noise power level, we use the same formula given in [START_REF] Georgiou | Low power wide area network analysis: Can LoRa scale[END_REF]: σ 2 n = -174 + NF + 10 log(B), where NF is the noise figure, usually set as 6 dB, B is the bandwidth occupied by the transmitted packet. The occupied transmission band is assumed to be 200 kHz (e.g. LoRaWAN). Thus, σ 2 n is -115 dBm. In terms of propagation path-loss model, Okumura-Hata model is applied. In urban area, the path-loss between one device with antenna height h d = 1.5 We consider a LPWAN network with target network packet loss rate 10%. With a given outage probability requirement, we calculate the minimum required BS spatial density (per km 2 ) with respect to the network load (per slot). The latter varies from 0.05 to 0.3. The numerical results are shown in Fig. 2 (with outage probability 10%) and Fig. 3 (with outage probability 1%).

From (15), we know that, if the background noise is neglected, there exists a linear relationship between network load (i.e. pλ m ) and minimum BS spatial density (i.e. λ b ) if target network packet loss rate and outage probability is given. This can be confirmed in Fig. 2 and Fig. 3 best Best BS attach with background noise, we unfortunately have no closedform expression between minimum required BS density and network load. From these two figures, the relationship between minimum required BS density and network load can be assumed to be linear and we can use curve fitting techniques to get a closed form expression for the case with noise with high accuracy. In addition, It is not surprised that with background noise, the minimum required BS density is more larger.

From these two figures, we observe that as network load increases, the increase of minimum BS spatial density in case of macro diversity is significantly smaller than that in case of best BS attach. One straightforward explanation is: to guarantee a certain level of outage probability, as the load increases, the network should be more dense and the distance between a considered device and its surrounding BS is closer. The chance that one packet is received by several surrounding BS is thus much higher in a dense network. Thus, macro diversity has more obvious effect.

Conclusion

In this paper, we study the performance of selection combining macro diversity in LPWAN with outage probability constraint. The outage probability is a performance indicator that better takes into account the real availability of the service in any part of the area. We compare the minimum required BS spatial density between best BS attach and SC macro diversity to confirm the interest of applying selection combining macro diversity. From numerical result, we observe that minimum required BS density and network load can be assumed to be linear if best BS attach method is applied and outage probability constraint is given. SC macro diversity is more resistant to such a constraint.

Fig. 1 .

 1 Fig. 1. Macro Reception Diversity Scheme Illustration

Fig. 2 .

 2 Fig. 2. Minimum BS spatial density with respect to network load. The target packet loss rate and outage probability threshold are both 10%.

Table 1 .

 1 Parameter values for numerical results BS with height h bs = 50 meters is[START_REF] Lagrange | Réseaux GSM: des principesa la norme[END_REF]: P ath Loss = -10 log 10 (K)+ γ10 log 10 (r) = 123.6 + 33.8 log(r), where log 10 (•) is 10-based logarithm operator, r is in unit of kilometer. Thus, the constant 10 log 10 (K) is -123.6 dB.

	Parameter Name	Parameter Notation Value
	Transmit power	Pt	14 dBm
	Noise figure	NF	6 dB
	Occupied bandwidth	B	200 kHz
	Target network packet loss rate P f, max	10%
	Outage probability	Poutage	1%, 10%
	Device height	h d	1.5 m
	BS antenna height	h bs	50 m
	Noise power	σ 2 n	-115 dBm
	meters and a

  . For slotted ALOHA based Minimum BS spatial density with respect to network load. The target packet loss rate is 10% and outage probability threshold is 1%.

		25	Best BS attach, slotted ALOHA									24.45 24.25
			Best BS attach, slotted ALOHA, no noise
	Mininum BS spatial density λb ) (per km 2	5 10 15 20	5.73 SC macro diversity, 10.19 1.06 1.69 4.85 Pure ALOHA with maximum interference 9.70 SC macro diversity, Pure ALOHA with maximum interference, no noise 2.12 2.57	14.55 14.88	3.52	3.18	19.65 19.40	4.51	4.24	5.53	5.30
		0	0.05	0.10	0.15 Network load pλm (per slot)	0.20	0.25
	Fig. 3.