Mattia Minelli

Maode Ma

Marceau Coupechoux

Jean-Marc Kelif

Marc Sigelle

Philippe Godlewski

Optimal Relay Placement in Cellular Networks

Keywords: Cellular network, relay, optimal placement, fluid model, simulated annealing, processor sharing

In this paper, we address the problem of optimally placing relay nodes in a cellular network with the aim of maximizing cell capacity. In order to accurately model interference, we use a dynamic framework, in which users arrive at random time instants and locations, download a file and leave the system. A fixed point equation is solved to account for the interactions between stations. We also propose an extension of a fluid model to relay based cellular networks. This allows us to obtain quick approximations of the Signal to Interference plus Noise Ratio (SINR) that are very close to 3GPP LTE-A guideline results in terms of SINR distribution. We then use these formulas to develop a dedicated Simulated Annealing (SA) algorithm, which adapts dynamically the temperature to energy variations and uses a combination of coarse and fine grids to accelerate the search for an optimized solution. Simulations results are provided for both in-band and out-of-band relays. They show how relays should be placed in a cell in order to increase the capacity in case of uniform and non-uniform traffic. The crucial impact of the backhaul link is analyzed for in-band relays. Insights are given on the influence of shadowing.

I. INTRODUCTION

Relaying is a promising feature of future cellular networks. The scenarios envisioned by the two standards IEEE 802.16j (for WiMAX networks) and 3GPP Release 10 and 11 (LTE-A) are the following: (a) coverage extension: relays should increase user experience in indoor or allow connection in shadowed zones; (b) group mobility: relays can aggregate the traffic related to a group of users within a train or a bus; (c) capacity boost: by deploying low-cost relay stations, a cellular operator can densify its network and increase its capacity. In this paper, we tackle the problem of optimal relay placement for capacity increase in an LTE-A-like cellular network.

A. Related Work

The relay placement problem arises in various contexts: wireless sensor networks (see e.g. [START_REF] Akkaya | Positioning of Base Stations in Wireless Sensor Networks[END_REF] and references therein), Wireless Local Area Networks (WLAN), WiMAX networks and cellular networks (like LTE-A); relays represent also a paradigm in information theory (see e.g. [START_REF] Chattopadhyay | Optimal Capacity Relay Node Placement in a Multi-hop Network on a Line[END_REF] and references therein). In this short literature review, we focus on the three Email addresses: mattia1@e.ntu.edu.sg (Mattia Minelli); emdma@ntu.edu.sg (Maode Ma); marceau.coupechoux@telecom-paristech.fr (Marceau Coupechoux); jeanmarc.kelif@orange.com (Jean-Marc Kelif); marc.sigelle@telecom-paristech.fr (Marc Sigelle); philippe.godlewski@godlewski@telecom-paristech.fr (Philippe Godlewski).

first cases characterized by a hierarchical structure made of Relay Nodes (RN) served by Base Stations (BS) or access points. The RNs location problem can be seen as a sub-case of the well studied facility location problem, which has been proved to be NP-hard [START_REF] Drezner | Facility Location: Applications and Theory[END_REF]. It deserves however some attention because the objective functions and the interdependence between transmitting stations through interference make it very specific.

Several papers in the literature focus on a single BS controlling several RNs [START_REF] So | Enhancing WLAN Capacity by Strategic Placement of Tetherless Relay Points[END_REF], [START_REF] Wang | Optimal Relay Location in Multi-Hop Cellular Systems[END_REF], [START_REF] Dong | Optimal Relay Location in OFDMA Based Cooperative Networks[END_REF], [START_REF] Chang | A Novel Relay Placement Mechanism for Capacity Enhancement in IEEE 802.16j WiMAX Networks[END_REF], [START_REF] Yang | Relay Station Placement for Cooperative Communications in WiMAX Networks[END_REF]. In [START_REF] So | Enhancing WLAN Capacity by Strategic Placement of Tetherless Relay Points[END_REF], authors consider the problem of optimal relay placement in a single WLAN cell. The problem consists in minimizing the average packet transmission time. The method is based on a discretization of the possible relay locations, a Lagrangian relaxation and an iterative algorithm. In [START_REF] Wang | Optimal Relay Location in Multi-Hop Cellular Systems[END_REF], RNs are constrained to be on a circle around the BS, so that authors come up with a single variable optimization problem. The same approach is taken in [START_REF] Dong | Optimal Relay Location in OFDMA Based Cooperative Networks[END_REF] except that cooperative strategies (Decode and Forward, DF and Amplify and Forward, AF) are assumed. In [START_REF] Chang | A Novel Relay Placement Mechanism for Capacity Enhancement in IEEE 802.16j WiMAX Networks[END_REF], a region controlled by a BS is divided into sub-regions characterized by their traffic requirements. The goal of the authors is to place relays such that capacity is maximized and user minimum bandwidth requirements are satisfied. Although useful, for example in rural areas, these approaches cannot be directly applied in a cellular network made of a dense network of BSs and RNs.

References [START_REF] Niyato | Relay-Centric Radio Resource Management and Network Planning in IEEE 802.16j Mobile Multihop Relay Networks[END_REF], [START_REF] Lu | Joint Base Station and Relay Station Placement for IEEE 802.16j Networks[END_REF], [START_REF] Yu | Planning Base Station and Relay Station Locations for IEEE 802.16j Network with Capacity Constraints[END_REF], [START_REF] Islam | Capacity-Optimal Relay and Base Station Placement in Wireless Networks[END_REF], [START_REF] Saleh | Analysis of the Impact of site Planning on the Performance of Relay Deployments[END_REF], [START_REF] Sambale | Cell Spectral Efficiency Optimization in Relay Enhanced Cells[END_REF], [START_REF] Wang | Approximation Algorithms for Cellular Networks Planning with Relay Nodes[END_REF] consider several BSs and several RNs per BS. A hierarchical optimization problem is formulated in [START_REF] Niyato | Relay-Centric Radio Resource Management and Network Planning in IEEE 802.16j Mobile Multihop Relay Networks[END_REF] for WiMAX networks: authors first focus on short term call admission control decisions (they use here the Markov Decision Process framework) and then, on the long term, on network planning (a binary integer linear problem, ILP, is solved with standard methods). Authors of [START_REF] Lu | Joint Base Station and Relay Station Placement for IEEE 802.16j Networks[END_REF] consider the joint deployment of BSs and RNs and try to maximize network capacity with a fixed budget. They formulate the problem as an ILP and propose a two stage RN and BS deployment algorithm to obtain sub-optimal strategies. A similar idea is used in [START_REF] Yu | Planning Base Station and Relay Station Locations for IEEE 802.16j Network with Capacity Constraints[END_REF]. Wang et al. [START_REF] Wang | Approximation Algorithms for Cellular Networks Planning with Relay Nodes[END_REF] try to minimize the installation cost for serving a given demand or to maximize the served demand for a fixed budget. Approximate algorithms are proposed and related approximation ratios are computed or bounded. Recently, [START_REF] Islam | Capacity-Optimal Relay and Base Station Placement in Wireless Networks[END_REF] has also tackled the joint problem of BSs and RNs placement, user allocation and transmit power setting. The goal is to maximize the network sum capacity while minimizing the installation cost. Sub-optimal solutions to the resulting mixed integer non linear program are obtained thanks to an iterative algorithm. These papers suffer from two main drawbacks. First, interference is never accurately modeled. It is out of the scope of [START_REF] Niyato | Relay-Centric Radio Resource Management and Network Planning in IEEE 802.16j Mobile Multihop Relay Networks[END_REF], [START_REF] Islam | Capacity-Optimal Relay and Base Station Placement in Wireless Networks[END_REF], [START_REF] Wang | Approximation Algorithms for Cellular Networks Planning with Relay Nodes[END_REF]; [START_REF] Lu | Joint Base Station and Relay Station Placement for IEEE 802.16j Networks[END_REF] considers only noise for the calculation of data rates; [START_REF] Yu | Planning Base Station and Relay Station Locations for IEEE 802.16j Network with Capacity Constraints[END_REF] tries to minimize the sum of path-losses, whereas data rate is related to the Signal to Interference plus Noise Ratio (SINR); [START_REF] Saleh | Analysis of the Impact of site Planning on the Performance of Relay Deployments[END_REF] assumes a single-interferer model. Second, user locations are often supposed to be known from the network planner (see e.g. [START_REF] Wang | Approximation Algorithms for Cellular Networks Planning with Relay Nodes[END_REF]). This assumption is realistic in a WiMAX network but cannot be considered for a cellular network.

Three recent papers take into account co-channel interference in their problem formulation. In [START_REF] Khakurel | Optimal Relay Placement for Coverage Extension in LTE-A Cellular Systems[END_REF], authors propose a fixed point algorithm for the RN placement problem. RNs are however constrained to be on a circle around their serving BS. [START_REF] Hamdi | Optimized Dual Relay Deployment for LTE-Advanced Cellular Systems[END_REF] relies on extensive system level simulations and constraints also RNs to be placed according to a predefined pattern. Sambale and Walke [START_REF] Sambale | Decode-and-Forward Relay Placement for Maximum Cell Spectral Efficiency[END_REF] propose a Simulated Annealing (SA) algorithm based on Monte Carlo simulations compliant with 3GPP LTE-A guidelines (although shadowing is ignored in the optimization). [START_REF] Guo | Relay Deployment in Cellular Networks: Planning and Optimization[END_REF] deals with RN positioning for network planning and optimization. Authors set up an analytic performance evaluation model, based on simple pathloss based SINR one-dimensional calculations with a single interferer. [START_REF] Sambale | Cell Spectral Efficiency Optimization in Relay Enhanced Cells[END_REF], which is an extension of [START_REF] Sambale | Decode-and-Forward Relay Placement for Maximum Cell Spectral Efficiency[END_REF], also copes with the problem of RN placement, taking into account physical and MAC layer in their model (shadowing is taken into account), and optimizing RNs placement in terms of cell spectral efficiency. However, this paper, as well as [START_REF] Khakurel | Optimal Relay Placement for Coverage Extension in LTE-A Cellular Systems[END_REF], [START_REF] Hamdi | Optimized Dual Relay Deployment for LTE-Advanced Cellular Systems[END_REF], [START_REF] Sambale | Decode-and-Forward Relay Placement for Maximum Cell Spectral Efficiency[END_REF], [START_REF] Guo | Relay Deployment in Cellular Networks: Planning and Optimization[END_REF] and other works, assumes that RNs and BSs transmit at every time instant (full buffer traffic). Hence, it does not model buffers loads and is equivalent to a static approach. This often results in misleading performance evaluation, as this model tends to overestimate interference.

B. Contributions

The contributions of this paper are the following:

• We introduce a "dynamic" framework for relays placement performance evaluation based on traffic analysis. This model overcomes the limits of widely used static models by considering users arriving in the system at random time instants and locations, downloading a file and leaving the system, rather than a fixed set of known user locations. RNs and BSs are modeled as M/G/1/PS queue and a fixed point iteration captures the interactions between transmitting nodes. Contrary to most of the literature, our framework takes into account non uniform traffic patterns. • We propose an extension of the fluid model, developed in [START_REF] Kelif | Downlink Fluid Model of CDMA Networks[END_REF] for interference modeling, to relay-based cellular networks in order to obtain quick calculations of the SINR. We thus take into account co-channel interference of the whole network in our optimization. Numerical results show that we are very close to 3GPP LTE-A guidelines results in terms of SINR distribution. • Based on these quick SINR calculations, we develop a dedicated SA algorithm for the relay placement problem.

Our algorithm enjoys some enhancements with respect to the standard SA. These enhancements have been already investigated in the field of image processing but have not been yet considered for wireless network optimization, to the best of authors' knowledge. More specifically, we propose a method to dynamically adapt temperature to energy variations and a combination of coarse and fine grids to accelerate the search for an optimal solution. The paper is organized as follows. In Section II, we present the system model, while in Section III, we deal with the performance of RNs placement, by formulating cell capacity. The SA algorithm is presented in Section IV and simulation results in Section V. Section VI concludes the paper.

II. SYSTEM MODEL

A. Network Topology and Serving Station Assignment

We consider a single frequency cellular network (all stations use the same frequency) consisting of tri-sectorized hexagonal cells. Every sector is controlled by a sector BS. In every hexagonal cell, n non-cooperative DF RNs [START_REF] Guo | Relay Deployment in Cellular Networks: Planning and Optimization[END_REF] are deployed. Every station (sector BS or RN) can be either active (i.e., transmitting) or idle (i.e., not transmitting), at any given time instant. We denote with P the transmit power of an active sector BS, and with P R the transmit power of an active RN. Let B and R be respectively the set of sector BSs and the set of RNs. We focus on capacity evaluation for the downlink.

The generic relay deployment is illustrated in Fig. 1. The deployment pattern is identical in all cells1 (the relative locations of the relays w.r.t. the cell center are constant across the cells). We label RNs of each cell with indices 1 • • • n, and define as type i relays those RNs labeled with i. The set of type i relays forms a regular stations pattern.

Each User Equipment (UE) is connected to its best server, i.e., the station which provides the highest signal power. We define S i as the region, of surface S i , where station i is the best server. Moreover, denoting with K c the set of stations belonging to a given hexagonal cell c, we define A c i∈Kc S i as the region, of surface A c , where users are served by a station of c.

B. Resource Organization

We consider in-band and out-of-band half-duplex relays. We assume a time division access between sector BSs and RNs illustrated in Fig. 2. Without loss of generality, the frame duration and the overall band are set to 1 time unit. When in-band relays are considered, a sector BS transmits data to the RNs it controls over the Backhaul Link (BL) during a time τ . Then, during 1τ and simultaneously, the BS transmits over the Direct Link (DL, BS-UE link) and RNs transmit over the Relay Link (RL, RN-UE link) to their respective attached UEs (Fig. 3). When out-of-band relays are considered, the BL is using another frequency band (e.g. over a microwave link) or a dedicated narrow beam, so that BL radio resources are not subject to the cellular network planning. In this case, τ = 0. These definitions of in-band and out-of-band relays are in accordance with the definitions of the 3GPP [START_REF][END_REF].

C. Propagation Model

Let consider a transmitting station j (RN on the RL, sector BS on the DL or BL) and a receiver u (UE on the DL or RL, RN on the BL) located in (r j , θ j ), where r j is the receiver-station distance and θ j is defined as the angle between the receiver-station direction and the station antenna boresight direction (θ j = 0 for omnidirectional antennas). We denote with superscripts L and N the propagation parameters referred to Line of Sight (LOS) and Non Line of Sight (NLOS) propagation respectively, and with subscripts B, R and D the propagation parameters referred to BL, RL and DL respectively. According to the 3GPP guidelines for relay performance evaluation [START_REF]Type-1 Relay Performance for Downlink[END_REF], the path-gain g l between a station and a receiver location can be written as follows: where l ∈ {B, R, D} depends on the link type, K L l and K N l are propagation constants, η L l and η N l are path-loss exponents, δ l (r j ) is a Bernoulli random variable (RV), which is equal to 1 when the considered link is LOS. The probability for a given link to be LOS depends on l and r j (for details see [START_REF]Type-1 Relay Performance for Downlink[END_REF]). Antenna pattern is given by A l (θ). RNs are equipped with omnidirectional antennas, hence A R (θ j ) = 1 ∀θ j . A sector BS uses the same antenna for BL and DL, so that A B = A D A. Shadowing is modeled by a log-normal RV X (u,j) l with standard deviation σ l . Note that the LOS or NLOS condition is supposed to be constant in time. Shadowing is assumed to change slowly in time, this notion will be detailed later in the paper.

g l (r j , θ j ) = δ l (r j )h L l (r j , θ j )+(1-δ l (r j ))h N l (r j , θ j ), (1)
= δ l (r j )A l (θ j ) K L l r η L l j +(1-δ l (r j ))A l (θ j ) K N l r η N l X (u,

D. Traffic Model

In this paper, we assume that each station is equivalent to an M/G/1/PS queue [START_REF] Kelly | Reversibility and Stochastic Networks[END_REF], [START_REF] Bonald | Wireless Downlink Data Channels: User Performance and Cell Dimensioning[END_REF]. This corresponds to a fair radio resource scheduling policy. Flow calls from users arrive according to a Poisson process of intensity λ(s)ds [f lows/s] in the location s of surface ds and flow calls sizes are i.i.d. with mean π [bit/f low]. We suppose for simplicity that the function λ(s) follows the same pattern in every cell. Now, let decompose λ(s) λφ(s), where the constant λ = Ac λ(s)ds/A c [f lows/s/m 2 ] is the average flow call intensity in A c , and φ(s) = λ(s)/ λ is named normalized flow intensity. The term φ(s) allows us to consider spatially nonuniform traffic and can be seen as a traffic profile. Moreover, we define ω(s) = λ(s)π [bit/s/m 2 ] as the traffic density in s. Finally, ω = λπ [bit/s/m 2 ] denotes the average traffic density in A c . The latter can be decomposed according to the surface controlled by each station: ω = n+3 i=1 ωi Si Ac , where ωi is the average traffic density on the surface S i controlled by device i: ωi π Si Si λ(s)ds. We denote with C(s) > 0 the user spectral efficiency (in bits/s/Hz) in location s, and with ρi the load of the M/G/1/PS associated to station i. If ρi > 1, station i is said to be saturated. According to the M/G/1/PS results, when ρi < 1, ρi is also the probability that station i is active (and thus interfere). By extension, if a station is saturated, it is always active. If α i is the Bernoulli RV, which is equal to 1 when station i is active and 0 otherwise, we have

P[α i = 1] = E[α i ] = min{ρ i , 1}.
By symmetry, all stations of the same type (sector BSs with the same antenna boresight direction or type i RNs) have the same load. At a given time instant t, the spectral efficiency in s, which depends on the SINR, is an explicit function of the α j , j = i, where i is the serving station, because of the interference term. We thus write C(s, α -i (t)), where α -i (t) is the vector of the activity variables of all stations in the network at time t, except i.

III. RELAYS PLACEMENT PERFORMANCE

Performance of RNs placement is measured with the corresponding cell capacity. This section is devoted to the assessment of the capacity of a cell made of three sectors and n RNs. Cell capacity (derived in III-C) depends on the BL capacity (III-A), on the spectral efficiency (III-D) and on the maximum traffic intensity stations can sustain, which is obtained by deriving stations loads (III-B). Subsection III-E introduces a method to speed up capacity evaluation, by approximating interference from far stations.

A. Backhaul Link Dimensioning

We derive here the proportion τ of radio resources dedicated to the BL, based on average RN throughputs, supposing that stations are not saturated. Let us consider a sector BS b and the set of p b RNs controlled by b (see Fig. (i,b) . Note that τ b is the load of the backhaul link between b and its corresponding RNs. In order to ensure a time synchronization between all sectors of a cell, we impose τ = max b {τ b }, so that in the best case:

i is not overloaded if C BL (i, b)τ i,b > ωi S i . Hence, a lower bound on τ b is given by: τ b > p b i=1 ωiSi C BL
τ max b p b i=1 ωi S i C BL (i, b) . (2)
Suppose that the flow call intensity is uniform, i.e., ω(s) = ω ∀s. In this case (2) simplifies to:

τ (ω) ω max b p b i=1 S i C BL (i, b) . (3)
In the following, we write τ (ω) to designate (2) or (3) in order to show the dependency of τ on the input traffic. The value of τ is set to zero when out-of-band relaying is adopted.

B. Stations Loads

The stations loads, which are related to their probability to be active, are coupled through the spectral efficiency function, which in turn is a function of the stations activities. We solve this problem using a fixed point iteration.

Lemma 1. The load of a station i is expressed by:

ρi = ω 1 -τ (ω) Si φ(s)E α 1 C(s, α -i ) α i = 1 ds. ( 4
)
Proof:

ρi

= lim

T →∞ 1 1 -τ (ω) T 0 Si ω(s) C(s, α -i (t)) 1l αi(t)=1 ds dt T 0 1l αi(t)=1 dt , (2)
= lim

T →∞ ω 1 -τ (ω) 1 T T 0 Si φ(s) C(s, α -i (t)) 1l αi(t)=1 ds dt 1 T T 0 1l αi(t)=1 dt , (3)
= ω 1 -τ (ω) E α Si φ(s) C(s, α -i ) 1l αi=1 ds P α (α i = 1) , (4)
= ω 1 -τ (ω) E α Si φ(s) C(s, α -i ) ds α i = 1 , (5)
= ω 1 -τ (ω) Si E α φ(s) C(s, α -i ) α i = 1 ds.
(1) comes from the fact that at a given instant t, the load generated in a small area ds in s is

ω(s)
C(s,α-i(t)) ds when the serving station i is active. Averaging is done on the time activity of station i. (2) comes from ω(s) = ωφ(s). (3) comes from the assumption that the process α(t) has a time-stationary limit α. (4) is by definition of the conditional expectation. ( 5) comes from the assumption that the shadowing changes slowly in time (and hence with respect to the realizations of α(t)) 2 . Hence, we assume that the expectation is taken over a period of time large enough for the steady-state of the M/G/PS/1 queues to be reached, but sufficiently small to keep the channel constant at each location. Recall moreover that the LOS and NLOS conditions are fixed in time. As a consequence, S i is fixed with respect to the realizations of α(t).

We denote ρ (ρ 1 , ..., ρn+3 ) the vector of loads corresponding to the n + 3 stations (n relays and three sector BSs) in the central hexagonal cell. For a given normalized flow intensity φ, let define the operator F (ρ, ω) = (F 1 (ρ, ω), ..., F n+3 (ρ, ω)) as follows:

F i (ρ, ω) = ω 1 -τ (ω) Si φ(s)E α 1 C(s, α -i ) α i = 1 ds, (5) where ∀i, P α [α i = 1] = E α [α i ] = min{ρ i , 1}. Let also define ωmax i ω 1-τ (ω) Si φ(s)ds C(s,1) . Theorem 1. If α → C(s, α) is a continuous mapping and 1/C(s, α) is non-decreasing in α, F : n+3 i=1 [0; ωmax i ] →
F i (ρ, ω) (1) = ω 1 -τ (ω) Si φ(s) lim N →∞ 1 N N h=1 1 C(s, α -i (h)) ds, (2)
≤ ω 1 -τ (ω) Si φ(s) lim N →∞ 1 N N h=1 1 C(s, 1)
ds,

(3) = ω 1 -τ (ω) Si φ(s) C(s, 1)
ds, A similar approach is considered in [START_REF] Rong | Performance Evaluation of Cellular Networks Offering TV Services[END_REF] (Corollary 1) for a cellular network without relays. We extend here the result to relay-based sectorized networks and we adopt a different mapping F that guarantees the existence of a fixed point (which is possibly outside the [0, 1) n+3 interval). As in [START_REF] Rong | Performance Evaluation of Cellular Networks Offering TV Services[END_REF], we cannot conclude on the uniqueness of the fixed point. Starting from a central isolated cell (i.e., without interference) in the fixed point iteration makes however sense in a context of increasing traffic.

The fixed point of ( 5) yields the loads ρi of all stations i = 1 • • • n + 3, for a given RN placement, a given normalized flow call intensity φ(s), and a given traffic density ω.

C. Cell Capacity

The capacity of the cell is defined as

C cell ωmax A c [bit/sec/Hz/cell], (6)
where ωmax is the maximum average traffic density that can be supported by a RN placement (without any station being saturated):

ωmax = max{ω ∈ R + : F (ρ, ω) = ρ and ρ ∈ [0, 1) n+3 }(7)
The value of ωmax can be found for any given normalized flow call intensity φ(s), by solving the fixed point of ( 5) for several values of ω and choosing the highest ω, for which no station is saturated, according to the desired accuracy. This can be done e.g. with a dichotomic search over ω.

D. Spectral Efficiency and SINR

We will now assume that the spectral efficiency in s is derived by means of a saturated Shannon formula:

C(s, α -i ) = min log 2 (1 + γ s (α -i )), C , (8)
where C is the maximum achievable spectral efficiency and γ s (α -i ) is the SINR in s. Note that this function fulfills the conditions of Theorem 1. Now, consider a UE u located in s and receiving from its serving station i a useful signal power P i (s) = Interference is exactly computed for near stations, while it is approximated for far stations by means of a fluid model of far stations network.

max j∈B∪R P j (s). The SINR γ s (α -i ) experienced by u is:

γ s (α -i ) = P i (s) I s (α -i ) + N th , (9)
where N th is the thermal noise power and I s (α -i ) = j∈B∪R, j =i α j P j (s).

E. Interference Computation with Fluid Model of Far Network Stations

Relays placements performance assessment through ( 6) and (8) involves intensive use of SINR evaluation: interference sum and SINR must be evaluated on the whole cell surface for each realization of α, at each iteration of the fixed point and for each proposed ω.

In this subsection, we propose a fast methodology for the computation of SINR, which should be simplified as much as possible in order to reduce the computation time of the optimization, while remaining accurate. Our approach decomposes the interference into two parts: the contribution of the near stations and that of the far stations. Far interference is approximated thanks to a fluid model [START_REF] Kelif | A Fluid Model for Performance Analysis in Cellular Networks[END_REF] adapted to relaybased cellular networks, while interference from near stations is exactly computed.

Let define the near BSs as the set B c of BSs in the central cell and in the first cells ring around the central cell, and far BSs as the set B f of BSs which are not near, so that B = B c ∪ B f and B c ∩ B f = ∅. Sets R c and R f are similarly defined for RNs (see Figure 4). Now, the overall interference I s (α -i ) can be decomposed into two contributions:

I s (α -i ) = I s,c (α -i ) + I s,f (α -i ), where I s,c (α -i ) = j∈Bc∪Rc, j =i α j P j (s) (10)
is the power received from near stations, while

I s,f (α -i ) = j∈B f ∪R f ,j =i α j P j (s) (11)
is the power received from far stations.

In our study, we compute I s,c (α -i ) explicitly, while I s,f (α -i ) is approximated with its average value over activity (cf. α) and shadowing (cf. X) variations:

I s,f (α -i ) ≈ E α,X [I s,f (α -i )] , (12)
= j∈B f E α,X [α j P j (s)] + k∈R f E α,X [α k P k (s)] .
In the above equation, expectations should be taken knowing α i = 1. However, we make the approximation that the activities of the far stations are independent on α i . It is reasonable to assume that a UE u in the central cell be always served by a near station, and that the propagation on all links between u and far stations be NLOS. This is valid considering commonly used shadowing standard deviations and LOS probability expressions (see e.g. [START_REF][END_REF], [START_REF]Type-1 Relay Performance for Downlink[END_REF]). Our simulations, based on the assumptions suggested in [START_REF]Type-1 Relay Performance for Downlink[END_REF], confirm these hypothesis, showing that less than 1% of users are connected to far stations, and LOS probability between u and a far station is close to zero. Hence, we assume I s,f (α -i ) to be composed by NLOS interferers, and express it as

I s,f (α -i ) ≈ j∈B f e aσ 2 D /2 E α α j P K N D r -η N D s,j A(θ) + k∈R f e aσ 2 R /2 E α α k P R K N R r -η N R s,k = e aσ 2 D /2 3 k=1 min{1, ρk } j∈B f,k P K N D r -η N D s,j A(θ) + e aσ 2 R /2 n h=1 min{1, ρh } k∈R f,h P R K N R r -η N R s,k , (13)
where r s,j is the distance between station j and s, B f,k is the set of far BSs of type k, R f,h is the set of RNs of type h and a = ln(10)/10. Note that parameters 13) are referred to NLOS propagation.

K N D , η N D , K N R and η N R in (
Sums j∈B f,k P K N D r -η N D s,j A(θ) and k∈R f,h P R K N R r -η N R s,k
can be approximated by adopting a fluid model [START_REF] Kelif | Downlink Fluid Model of CDMA Networks[END_REF] for the far stations networks

B f,k , k ∈ {1 • • • 3} and R f,h , h ∈ {1 • • • n}.
The fluid model is a powerful tool to simplify interference computation in hexagonal [START_REF] Kelif | A Fluid Model for Performance Analysis in Cellular Networks[END_REF] and dense Poisson [START_REF] Combes | A Justification of the Fluid Network Model Using Stochastic Geometry[END_REF] cellular networks. It is based on approximating a discrete set of network stations, lying on a given region, with a continuum of stations. The density (measured in [stations/m 2 ]) of the continuum is set to be equal to the density of discrete stations in the original network. Interference sum is then approximated by integrating received power from the continuum of stations, over the considered area. The main advantage of this approach is that it allows us to obtain approximate closed-form interference expressions, which solely depend on the distance between s and the network center. We refer the reader to [START_REF] Kelif | A Fluid Model for Performance Analysis in Cellular Networks[END_REF] for a detailed explanation and validation through Monte Carlo simulations.

Let now focus on

j∈B f,k P K N D r -η N D s,j A(θ), k ∈ {1 • • • 3}.
We adopt the fluid model, substituting the hexagonal network of type-k far BSs with a continuum of BSs of the same type, lying on a ring centered at s. Following the approach presented in [START_REF] Kelif | Fluid Model of the Outage Probability in Sectored Wireless Networks[END_REF], interference sum can be approximated as

j∈B f,k P K N D r -η N D s,j A(θ) ≈ ρ BS ∞ Rc-r0 2π 0 P K N D r -η N D A(θ) rdrdθ, = ρ BS b P K N D η N D -2 (R c -r 0 ) 2-η N D , (14)
where ρ BS is the BS sites density, b = 2π 0 A(θ)dθ, r 0 is the distance between the location s and the cell center and R c is the distance between the closest far BS and the cell center.

Terms

k∈R f,h P R K N R r -η N R s,k , h ∈ {1 • •
• n} can also be approximated adopting a fluid network model, obtaining

k∈R f,h P R K N R r -η N R s,k ≈ ρ BS ∞ Rc-r h 2π 0 P R K N R r -η N R rdrdθ, = 2πρ BS P R K N R η N R -2 (R c -r h ) 2-η N R , (15)
where r h is the distance between s and the RN of type h in the central cell. Using ( 14) and ( 15), I s (α -i ) is finally approximated as

I s (α -i ) ≈ j∈Bc∪Rc, j =i α j P j (s) (16)
+e aσ 2 D /2 3 k=1 min{1, ρk }ρ BS b P K N D η N D -2 (R c -r 0 ) 2-η N D +e aσ 2 R /2 n h=1 min{1, ρh }2πρ BS P R K N R η N R -2 (R c -r h ) 2-η N R .
The use of a fluid model allows to approximate the sum of interference from all far stations, at each considered cell point, by only computing n + 1 closed-form formulae, while the hexagonal model requires to calculate the sum of (n + 3) N R i=N f 6i far stations received powers, where N R denotes the number of cells rings around the central cell and N f the first ring of far cells. This considerably reduces computational burden, considering that interference calculation is inserted in multiple nested computation cycles, as mentioned above.

IV. OPTIMAL PLACEMENT

In this section, we study the optimization of relays placement, by means of a dedicated Metropolis-Hastings Simulated Annealing (SA) optimization algorithm, with the target of maximizing cell capacity C cell , defined in Section III (see [START_REF] Dong | Optimal Relay Location in OFDMA Based Cooperative Networks[END_REF]). We first detail the configuration space, i.e., the set of variables to be sought, and then describe the used SA algorithm, with several enhancements which have been implemented.

In this work, a configuration is given by the positions of all relays in a cell, which we assume to lay on an hexagonal grid, due to the symmetry of the problem (see Fig. IV-1). As the RN location problem includes the capacitated facility location problem as a special case, it is NP-hard [START_REF] Drezner | Facility Location: Applications and Theory[END_REF]. For a typical grid of size N = 1024 measurements points per cell and n = 6 relays, the cardinality of the configuration space Ω is:

|Ω| = C N n = C 1024 6
> (1000) 6 /6! > 10 15 . Hence, exhaustive search is infeasible in practice and we have to address other optimization techniques. Simulated Annealing (SA) is a wellknown stochastic technique for solving such large combinatorial optimization problems endowed with fairly non-trivial energy landscapes. It originates to [START_REF] Metropolis | Equation of State Calculations by Fast Computing Machines[END_REF] but was rediscovered later [START_REF] Laarhoven | Simulated Annealing: Theory and Applications[END_REF], and with great success in network optimization up to now [START_REF] Keung | The Base Station Placement for Delay-Constrained Information Coverage in Mobile Wireless Networks[END_REF], [START_REF] Chen | Gibbsian Method for the Self-Optimization of Cellular Networks[END_REF]. It leads to efficient optimization if its parameters (especially the temperature schedule) are well setup [START_REF] Robini | Simulated Annealing, Acceleration Techniques, and Image Restoration[END_REF]. Recall that Ω is a finite configuration space and that we consider a cost energy function U (x) : Ω → R to be minimized. In this work, the energy of a candidate configuration x is the inverse of its related cell capacity (see ( 6)):

U (x) = -C cell (x). (17)
The principle of SA lays in assigning the following exponential probability to any configuration:

P(x) = e -U (x) Z ∀x ∈ Ω (with Z = x∈Ω e -U (x) ). (18)
A minimizer of U (.) possesses thus maximal probability P(.) and can be found as follows (using Metropolis-Hastings variant [START_REF] Hastings | Monte Carlo Sampling Methods using Markov Chains and their Applications[END_REF]):

• Initialization: assign an arbitrary initial configuration x 0 ∈ Ω.

• At step m ≥ 0: let x = x m be the current configuration. Apply the following procedure: pick up a candidate configuration x ′ ∈ Ω according to a user-specified proposal law r(x → x ′ ) and compute then the following acceptance rate:

Ξ(x → x ′ ) = min 1, P(x ′ ) P(x) 1 Tm . r(x ′ → x) r(x → x ′ ) , = min 1, e -U (x ′ )-U (x) Tm . r(x ′ → x) r(x → x ′ ) . (19)
Assign x m+1 = x ′ with probability p = Ξ(x → x ′ ) (and

x m+1 = x with probability 1 -p).
Here, T m is a positive temperature parameter depending on step m required to slowly decrease to 0 as m → +∞ (and rigorously to satisfy T m ≥ T0 1+log(m+1) ). Usually a geometric law is adopted T m = T 0 . β m with β < 1 but close to 1. Hereafter an adaptive temperature schedule is investigated. Notice also that when the proposal law is uniform or more generally symmetric, i.e., r(x ′ → x) = r(x → x ′ ), the acceptance rate boils down to the usual Metropolis form:

Ξ(x → x ′ ) = min 1, e -U (x ′ )-U (x) Tm . ( 20
)
This generic method enjoys a number of extensions and variants that are implemented here: 1) "Restricted Image Spaces" [START_REF] Robini | Simulated Annealing, Acceleration Techniques, and Image Restoration[END_REF]: It is indeed preferable to draw at each step a configuration x ′ which is close to x, for instance by varying one variable only (say, the position ξ i of the relay of type i) and in a restricted range around its current ξ i ξ' i Fig. 5.

Typical hexagonal cell with a hexagonal grid of possible relay locations. The current location is ξ i and the acceptance law is Gaussian along a random direction. If the candidate ξ ′ i falls outside the hexagon, it is "periodized" inside the cell.

value. This is known to increase the global algorithm search speed [START_REF] Robini | Simulated Annealing, Acceleration Techniques, and Image Restoration[END_REF]. To this aim we adopt a continuous framework, with symmetric gaussian proposal probability distribution function (PDF) given by:

r(ξ i , ξ ′ i ) = 1 2π 1 ζ 2 N ( ξ ′ i -ξ i ; ζ 2 ). (21)
Here, ζ controls the average distance between sites ξ i and ξ ′ i (Fig. IV-1). Finally the closest discrete grid site to ξ ′ i is selected. This is justified (mostly at fine resolution) for reasonable values of ζ w.r.t. the lattice step: first, there is low probability to reject proposed grid node, i.e., to find an already busy site. Then, the overall periodicity of RN lattice is injected by periodizing the proposed site itself into the whole lattice (Fig. IV-1). Thus, with this approximation the proposal laws cancels down in Hastings-Metropolis, and we obtain the usual Metropolis acceptance ratio.

2) Adaptive temperature scheme: In this work, we try to automatically (and adaptively) estimate both initial temperature T 0 and the cooling coefficient β. This has been previously done in the literature [START_REF] Robini | Simulated Annealing, Acceleration Techniques, and Image Restoration[END_REF] provided that one knows the energy landscape (which is not the case here). This is done in two steps:

• Initial temperature setting: it is adjusted by imposing that the initial average acceptance ratio

Ξ 0 = <Ξ(x → x ′ )> ∈ [a, b] (in practice, we select [a, b] = [0.5, 0.8]).
To do this, T 0 is first kept constant during some number of steps (depending on the size of Ω). Then, a dichotomic-type update process is applied:

T 0 ← β T 0 if Ξ 0 > b and T 0 ← 0.5 β T 0 if Ξ 0 < a.
• Cooling schedule: after the previous initialization phase, a fixed number of simulated annealing steps M l and a final temperature T M l [START_REF] Robini | Simulated Annealing, Acceleration Techniques, and Image Restoration[END_REF] are assigned at each scale l = 1, 2 (see below), with related temperature schedule:

T m l = T 0 l × (T M l /T 0 l ) m/M l . ( 22
)
3) Multiscale implementation: Multiscale algorithms have been employed long ago in various branches of Applied Mathematics. In image processing several multiscale strategies have been outlined in [START_REF] Graffigne | Hierarchical Markov Random Field Analysis Models Applied to Image Analysis: a Review[END_REF], and fast algorithms for movement detection/segmentation of video sequences can be found in [START_REF] Pérez | Restriction of a Markov Random Field on a Graph and Multiresolution Statistical Image Modeling[END_REF]. For instance a hand movement can be decomposed first in a large range displacement (translation-rotation) of the hand, followed by detection of fingers moves at a finer scale. Similarly here we use a two-step setup with a first optimization on a coarse grid, then refinement of this solution at a finer level. Both steps employ SA with their own adapted parameters (initial temperature, temperature step and spatial grid resolution). The advantage of such a multiscale strategy seems twofold here. First, it allows, at lower level, to avoid spurious transitions between (high level) symmetry-invariant configurations, which can arise due to the complexity of the global energy landscape. Then, it enables to check the robustness at both coarse and fine level of the placement solution to shadowing, LOS and NLOS conditions.

V. RESULTS

A. Simulation Assumptions 1) Spectral Efficiency Computation:

Simulation assumptions related to network layout and propagation (see ( 1)) follow the setup described in [22, Appendix: Simulation Assumptions] (case 3) for LTE-Advanced. In particular, the sector BS and RN transmit powers are resp. 46 dBm and 30 dBm and the cell range is 1 km. The shadowing standard deviations are 10, 8 and 6 dB on the RL, DL, and BL resp. for NLOS propagation. Shadowing RVs between cell sites and between sector BSs are correlated (see [START_REF] Saleh | Analysis of the Impact of site Planning on the Performance of Relay Deployments[END_REF]). The correlation distance is 50 m (we use the Cholesky decomposition as in [START_REF] Klingenbrunn | Modelling Cross-Correlated Sahdowing in Network SImulations[END_REF]).

The SINR is computed on a regular hexagonal grid of Measurement Points (MP) in the cell (with 25 m spacing), in order to adapt to the problem geometry. For a given realization of the propagation random variables, a given value of ω, a given vector of station loads ρ and a given realization of the activity vector α, the SINR is computed according to [START_REF] Niyato | Relay-Centric Radio Resource Management and Network Planning in IEEE 802.16j Mobile Multihop Relay Networks[END_REF].

Instead of computing spectral efficiency using (8), we rely on the Modulation and Coding Schemes (MCS) indicated in [START_REF][END_REF] in order to have more realistic results. According to [START_REF][END_REF], below a certain SINR threshold, a MP is in outage, i.e., this position cannot be served by any station. In this case, the contribution of such MP to stations loads is zero. Any placement with an outage probability greater than 1% is rejected by the optimization process (see [START_REF] Fischione | Power and Rate Control with Outage Constraints in CDMA Wireless Networks[END_REF], [START_REF] Alsawah | Base-Station and Subcarrier Assignment in Two-Cell OFDMA Downlink under QoS Fairness[END_REF]).

2) Fixed Point Iteration: The expectation of the spectral efficiency is computed on every MP over hundred realizations of α and a new vector ρ is deduced from (4). We iteratively solve the fixed point equation F (ρ, ω) = ρ (about six iterations ensures convergence in our context). If ∃i, s.t. ρi ≥ 1, this means that ω > ωmax , otherwise ω < ωmax , see [START_REF] Chang | A Novel Relay Placement Mechanism for Capacity Enhancement in IEEE 802.16j WiMAX Networks[END_REF]. The values of ω to be used in the fixed point iteration are selected based on a dychotomous update process, which ensures a maximum error for C cell equal to ±2.3 × 10 -3 [bit/s] (this corresponds to ±23 [kbit/(s × cell)] for a system bandwidth of 10 MHz).

3) Simulated Annealing: Optimization of RN positioning is performed by evaluating the energy of a number of candidate RN placements during the execution of the SA. Each RN Fig. 6. SINR distribution, our SINR approximation vs. 3GPP results [START_REF]Type-1 Relay Performance for Downlink[END_REF] (static case, n = 0 or 3 relays).

is allowed to be located on a hexagonal grid spanning the whole cell. Two RNs are not allowed to be located on the same spot, or in the cell center. During the first optimization phase, RNs can be located on a coarse hexagonal grid (see Section IV), where the spacing between 2 neighboring MPs is 200 m. During the second phase, the spacing becomes 50 m. In this phase, candidate RN locations cannot be farther than 300 m from the location found at the end of the first phase. The number of states analyzed during each phase varies according to n. For example, if n = 3, the first phase of the SA is composed of M 1 = 30 temperatures steps.

During each temperature step, the energy of 250 candidate placements are evaluated. The second phase of the SA is characterized by M 2 = 20 and 150 placements for each step.

The SA is stopped if the acceptance rate Ξ is zero for two consecutive temperature steps. In this case, we assume that the algorithm has already reached a 'stable' solution, and we elect as final placement the one with the lowest energy among those analyzed up to the algorithm stop. Note that the implemented SA provides optimized solutions and not necessarily optimal solutions that can be obtained theoretically after an infinite number of iterations.

B. Model Validation

We first show the accuracy of SINR approximation introduced in Section III-E. Fig. 6 plots the Cumulative Distribution Function (CDF) of the SINR for the static case (α(t) = 1, ∀t) with n = 3 relays and without relays. We have compared the curves obtained in [START_REF]Type-1 Relay Performance for Downlink[END_REF] with those derived by means of our SINR approximation under the same assumptions. The results show that there is a good match between our fluid model approach and the results obtained by 3GPP. This can be explained by the shadowing model (the standard deviation is relatively low and RVs are correlated) and the resulting low probability for a UE to be attached to a far station.

C. Simulation Results

Unless specified, simulation assumptions are taken from [START_REF]Type-1 Relay Performance for Downlink[END_REF] (case 3). In Fig. 7, we first show the influence of the In-band (no shadowing)

Out-of-band (no shadowing)

In-band (with shadowing)

Out-of-band (with shadowing)

Out-of-band

In-band number of RNs and of their transmit power on their distance to the cell center. Let us first consider out-of-band RNs. The absence of backhaul constraints allows the RNs to cover the most interfered regions of the cell, i.e., the cell edge.

Increasing n tends to decrease the average distance to the cell center because of a repulsion effect: RNs interfere more and tends to move away from each other. While a set of RNs is remaining close to cell edge, another set moves closer to the BS site so that several rings of RNs may appear. Increasing P R induces a similar but smaller repulsion effect for n ≥ 5. When n is small, P R has however less influence because the inter-RN distances are bigger and thus their mutual influence is lower. The in-band case is characterized by a trade-off between the advantage of covering the cell edge and the price to be paid on the backhaul in terms of capacity. When P R is increased, a RN controls a bigger region (especially when n is small) and thus require a higher data rate on the backhaul. Consequently, RNs have to be placed closer to the sector BSs. In Fig. 8, we can observe that the shadowing has a negligible effect on the location of out-of-band RNs. The best server mechanism assumed on the DL and the RL indeed compensates the effect of the channel variations. The impact is however decisive for in-band RNs. The best server policy is indeed not assumed on the BL 3 . This means that the BL is possibly highly interfered and its capacity is degraded. This explains that RNs have to move closer to the sector BS in order to benefit from higher MCSs and LOS propagation.

Fig. 9 and 10 show some examples of optimized RN placements. In the out-of-band case, we see how RNs are preferably placed on the cell edge. With one or two RNs (not shown on the figure), the RNs cover corners of the hexagon. When the number of RNs increases, e.g. with 6 RNs, we see how three of them move closer to the cell center. If n still increases a second ring of RNs around the BS site appears. It is also noticeable from these figures that the RN placement in a given cell is coherent with the RN placement in neighboring cells. With four RNs for example, RNs attached to different cell sites form a regular pattern around the edge of the central cell. In the in-band case, RNs are much closer to the BS site because of the BL influence. The boresight direction of the sector BSs, which has almost no influence on the out-of-band case, plays now an important role. On the one hand, RNs tends to be in the boresight direction in order to benefit from a better backhaul. On the other hand, most interfered regions lies on the frontier between two sectors. For two RNs for example, the first effect is preponderant.

Fig. 11 shows the cell capacity as a function of n. In the outof-band case, increasing n or P R leads to a capacity increase. When there is no backhaul constraint, relaying is indeed equivalent to a classical network densification, which clearly increases the network capacity [START_REF] Kelif | Cell Breathing, Sectorization and Densification in Cellular Networks[END_REF]. Higher is P R , higher is the offload of the sector BSs towards RNs. As they control bigger regions, sector BSs are indeed the limiting stations: their load reaches 1 well before the RN loads do. According to [START_REF] Chang | A Novel Relay Placement Mechanism for Capacity Enhancement in IEEE 802.16j WiMAX Networks[END_REF], this limits in turn the cell capacity. Increasing P R consequently balances the traffic among RNs and BSs, which results in a higher capacity. For in-band RNs, there are two contradictory effects: increasing P R offloads the sector BSs but increases also the proportion of resources dedicated to the BL (RNs control bigger regions). Numerical results show that the first effect is slightly preponderant. The increase of capacity observed for out-of-band as well as in-band RNs contradicts the effect observed in [START_REF] Guo | Relay Deployment in Cellular Networks: Planning and Optimization[END_REF], which suggests a capacity decrease after 3 to 7 RNs (depending on P R ). This can be explained by the fact that [START_REF] Guo | Relay Deployment in Cellular Networks: Planning and Optimization[END_REF] assumes a full buffer traffic model (stations are always active) and supposes that RNs are placed on a circle around the BS. Fig. 12 confirms the influence of shadowing: it has a relative small impact on out-of-band RNs but greatly degrades the in-band RNs performance.

These results clearly show the advantage of deploying outof-band relays, especially in terms of cell capacity. The benefit of in-band relays in a scenario similar to the one defined by the 3GPP is less obvious. First, the increase of capacity is small. Second, RNs are placed close to the BS because of the backhaul constraints, so that the capacity increase benefits mainly to UEs having already good radio conditions without relays. The deployment of in-band RNs can thus be interesting only if the BL benefits from a much higher capacity either because of a good radio propagation (as in [START_REF] Combes | Self-Organizing Relays: Dimensioning, Self-Optimization, and Learning[END_REF]) or because RN cell selection is adopted (as in [START_REF] Saleh | Analysis of the Impact of site Planning on the Performance of Relay Deployments[END_REF]).

Finally, Fig. 13 shows an example of non-uniform traffic pattern inside the cell. The normalized traffic intensity φ(s) follows a bivariate normal distribution centered in (0.3, 0.4) km and with standard deviation 0.3 km. Iso-flow intensity levels are shown with concentric circles. As expected, the relays are placed closed to the hot spot. In the out-ofband case, the relays are also close to the cell border in order to cover the most interfered region. On the contrary, in-band relays are between the BS site and the hot spot in order to benefit from a good BL. Now the question arises of the accuracy of the solutions obtained by SA. There are few results on this problem in the literature and they are generally based on some consideration over the energy surface. For example, the work of Catoni (see e.g. [START_REF] Catoni | Sharp Large Deviations Estimates for Simulated Annealing Algorithms[END_REF]) is based on large deviations theory, and deals with the evaluation (in a probabilistic framework) of the distance of the final attained configuration to the optimal one. However, the performed analysis seems not to be applicable in our case, since we have a very weak knowledge of the energy landscape of this complex problem.

VI. CONCLUSION

In this paper, we have addressed the problem of the optimal placement of relays in a cellular network with the aim of increasing the downlink cell capacity. Compared to other works on this subject, traffic and SINR computation models are more realistic: we have set up a dynamic traffic model, where each station is equivalent to a M/G/1/PS queue. The interaction between stations is captured by a fixed point equation. The system is stable if none of the stations is overloaded. This approach accurately models the activity of the stations and the relative weight of most interfered regions. Uniform and non-uniform traffic can be analyzed. The optimization of the placement is done using a dedicated Simulated Annealing algorithm. In order to speed up the search for an optimized solution, we have developed an extension of a fluid model In-band (no shadowing)

Out-of-band (no shadowing)

In-band (with shadowing)

Out-of-band (with shadowing)

In-band

Out-of-band Fig. 12. Cell capacity vs. n (in-band and out-of-band, P R = 30 dBm, with and without shadowing).

Hot spot center

In-band Out-of-band Hot spot center Fig. 13. Optimized out-of-band and in-band relay locations with a nonuniform traffic pattern (P R = 30 dBm, n = 3, the BS site is in the center of the hexagons, RNs are filled circles).

to relay based cellular networks. This approach yields very good approximations of the SINR at every location. Our SA algorithm adapts dynamically the temperature to the energy variations and uses a combination of coarse and fine grids. Simulation results shows that out-of-band relays are preferably placed on the cell edge and are arranged in rings around the BS, when their number increases. In-band relays suffer from the poor quality of the backhaul link, especially in presence of shadowing and tends to be much closer to the BS. In both cases, cell capacity increases with the number of relays. The benefit is however small with in-band relays.

Fig. 1 .

 1 Fig. 1. Example of relay deployment with n = 4 relays. The cell range is R. Sites are tri-sectorized and boresight directions of sector BSs are shown with arrows.

Fig. 2 .

 2 Fig. 2. Frame structure on the downlink for (a) in-band relays and (b) out-of-band relays when BS sector controls p ≤ n relays.

Fig. 3 .

 3 Fig. 3. Relay Link (RL), Backhaul Link (BL), and Direct Link (DL).

  2). Denote τ i,b the proportion of resources needed by b to serve RN i and τ b = p b i=1 τ i,b . The average throughput of RN i is given by Si ω(s)ds = ωi S i . Let C BL (i, b) be the average throughput on the BL between b and i. The link between b and RN

]

  has at least one fixed point. Proof: We have for any ρ ∈ n+3 i=1 [0; ωmax i ] and any station i, 0 ≤ F i (ρ, ω) and:

,

  where α -i (h) indicates the h-th realization of RV α -i . (1) comes from the definition of average. (2) results from the fact that, for each realization h of α -i , we have 1/C(s, α -i (h)) ≤ 1/C(s, 1). (3) follows from the consideration that all terms (1/C(s, 1)) in the summation have the same value. Finally, (4) is by definition of ωmax i . Now, using the Brouwer's fixed point theorem, we conclude the proof.

Fig. 4 .

 4 Fig. 4.Interference is exactly computed for near stations, while it is approximated for far stations by means of a fluid model of far stations network.

nFig. 7 .

 7 Fig. 7. Average RN to cell center distance vs. n, the number of RNs per cell (in-band and out-of-band, P R = 26, 30 and 46 dBm).

n

  Average RN-to-cell-center distance[km]

Fig. 8 .

 8 Fig. 8. Average RN to cell center distance vs. n (in-band and out-of-band, P R = 46 dBm, with and without shadowing).

Fig. 9 .Fig. 10 .

 910 Fig. 9. Optimized out-of-band relay locations (P R = 30 dBm, n = 2, 4 and 6, RNs are empty circles, BS sites are filled circles).

Fig. 11 .

 11 Fig. 11. Cell capacity vs. n (in-band and out-of-band, P R = 26, 30 and 46 dBm).

n

  Ccell [bit/s/H z/cel l ]

Mattia

  Minelli received his Master degree in telecommunications engineering from Politecnico di Milano in Italy, in 2009. He is currently participating in a Joint PhD Program between the School of Electrical and Electronic Engineering at Nanyang Technological University in Singapore and the Computer and Network Science Department of Telecom ParisTech in Paris. His research interests include wireless networks, propagation and radio resource management issues, and relaying in OFDMA networks. Ma Maode received his PhD degree in computer science from Hong Kong University of Science and Technology in 1999. Now, Dr. Ma is an Associate Professor in the School of Electrical and Electronic Engineering at Nanyang Technological University in Singapore. He has extensive research interests including wireless networking and network security. Dr. Ma has more than 250 international academic publications including over 100 journal papers and more than 130 conference papers. He currently serves as the Editor-in-Chief of International Journal of Electronic Transport. He also serves as an Associate Editor for other five international academic journals. Dr. Ma is an IET Fellow and a senior member of IEEE Communication Society and IEEE Education Society. He is the vice Chair of the IEEE Education Society, Singapore Chapter. He is also an IEEE Communication Society Distinguished Lecturer. Marceau Coupechoux has been working as an Associate Professor at Telecom ParisTech since 2005. He obtained his Masters' degree from Telecom ParisTech in 1999 and from University of Stuttgart, Germany, in 2000, and his Ph.D. from Institut Eurecom, Sophia-Antipolis, France, in 2004. From 2000 to 2005, he was with Alcatel-Lucent. He was a Visiting Scientist at the Indian Institute of Science, Bangalore, India, during 2011-2012. Currently, at the Computer and Network Science department of Telecom ParisTech, he is working on cellular networks, wireless networks, ad hoc networks, cognitive networks, focusing mainly on radio resource management and performance evaluation. Jean-Marc Kelif received the Licence of Fundamental Physics Sciences from the Louis Pasteur University of Strasbourg (France), and the Engineer degree in Materials and Solid State Physics Sciences from the University of Paris XIII in 1984. He worked as engineer in materials, and specialized in telecommunications. Since 1993, he has been with Orange Labs in Issy-les-Moulineaux, France. His current research interests include the performance evaluation, modeling, dimensioning and optimization of wireless telecommunication networks. Marc Sigelle graduated as an Engineer from Ecole Polytechnique in 1975 and from Telecom ParisTech in 1977. He obtained a PhD from Telecom ParisTech in 1993. He worked first at Centre National d'Etudes des Tlcommunications in physics and computer algorithms. Since 1989 he has been working at Telecom ParisTech in image processing, his main topics of interests being the restoration and segmentation of images using stochastic/statistical models and methods (Markov Random Fields, Bayesian networks, Graph Cuts in relationships with Statistical Physics). In the past few years he has been involved in the transfer of this knowledge to the modeling and optimization (namely, stochastic) of wireless networks. Philippe Godlewski received the Engineers degree from Telecom Paris (formerly ENST) in 1973 and the PhD in 1976 from University Paris 6. He is Professor at Telecom ParisTech in the Computer and Network Science department. His fields of interest include Cellular Networks, Air Interface Protocols, Multiple Access Techniques, Error Correcting Codes, Communication and Information Theory.

The possibility to change the deployment pattern from cell to cell could improve performance. This aspect is left for further studies.

Measurements for outdoor static scenarios reported in[START_REF] Seetharam | A Markov Chain Model for Coarse Timescale Channel Variation in an 802.16e Wireless Network[END_REF] justify our hypothesis, explaining this result with the fact that most of physical elements causing shadowing outdoor are fixed.

Reference[START_REF] Saleh | Analysis of the Impact of site Planning on the Performance of Relay Deployments[END_REF] proposes to attach a RN to the best BS, it is however not an option offered by the standard so far.

ACKNOWLEDGMENT

The authors would like to thank the anonymous reviewers for their comments and suggestions that helped improve the quality of this paper. M. Minelli is supported by the A-STAR Singapore International Graduate Award. This work is also supported by the french ANR project NETLEARN (ANR-13-INFR-004).