Olivier Roussel

Michel Taïx

Timothy Bretl

Motion Planning for a Deformable Linear Object

We present a motion planning algorithm for quasi-static Kirchhoff elastic rod in complex environments. By considering the set of quasi-static deformations which forms a finite dimensional manifold that can be parameterized by a single chart, the configuration space formulation extends nicely to this deformation space. As this parameterization is computationally expensive, our algorithm takes advantage of its smooth structure to approximate fast its neighborhood. Our approach enables motion planning for physically realistic deformable models for rods and our algorithm gives a significant performance improvement over classical motion planning approaches. We demonstrate the effectiveness of our approach on various academic and industrial scenarios, including difficulties such as narrow passages and high number of primitives.

I. INTRODUCTION

Motion planning in Robotics has been extensively studied from last decades. But most of the work focused on rigid bodies and relatively little work has been given to extension of the problem to a deformable robot. Recent applications of robotics algorithmic in various fields such as virtual prototyping have brought new motivations in this direction. It is impossible to correctly modelize industrial case (used in automotive and aeronautics industry) whitout deformable parts. These parts are not handled yet due to the lack of efficient motion planning algorithms and this is becoming a major limitation.

In this paper we consider only the case of "Deformable Linear Objects", (DLOs), robot in a rigid environment. DLOs are characterized by having one dimension much more greater than the two others (cable, hose, pipe,...).

For a rigid body, the configuration is defined by the finite number of parameters (here 6) that defines the frame, whereas for a deformable robot the number of shapes given a fixed frame is infinite. As deformation can be seen as adding new degrees of freedom to the system, the extension to current configuration space based formulation might seem straightforward. However, working with an inefficient parameterization of a deformation would lead to a high number of dependent degrees of freedom. This involves a high-dimensional configuration space although the set of deformations describes a much lower-dimensional manifold. Also, the deformation model that would describe how to get from the configuration space to the workspace may be hard to compute.

Flexible rods mechanics have been extensively studied, especially in the case of elastic deformations [START_REF] Stuart S Antman | Kirchhoff's problem for nonlinearly elastic rods(for axial extension and cross section shear)[END_REF] [START_REF] Simo | A finite strain beam formulation. the three-dimensional dynamic problem. part i[END_REF] [10] [START_REF] Bergou | Discrete Elastic Rods[END_REF]. These works offer relatively accurate dynamic models of an elastic rod, but the computational cost induced by the use of numerical methods such as finite elements is too high to include them directly into a motion planning framework. For instance, in [START_REF] Rodríguez | Planning motion in completely deformable environments[END_REF], Rodriguez et al. coupled a deformable dynamics simulator with a kynodynamic motion planning algorithm. However, the use of fully deformable environments prevents the robot to be stuck in local minima and bypasses the local control problem.

As in our context of assembling and disassembling studies the goal is to find a geometric path for the rod, the use of quasi-static models seems a reasonable assumption. In this direction, the work from Lamiraux and Kavraki [START_REF] Lamiraux | Planning paths for elastic objects under manipulation constraints[END_REF] investigated for the manipulation planning problem of deformable objects based on the computation of their equilibrium configurations. This computation relies on the numerical minimization of the total elastic energy for given gripper placements. Some efforts have been done by Moll and Kavraki [START_REF] Moll | Path planning for deformable linear objects[END_REF] in the special case of inextensible DLOs to reduce the cost of the minimization by the use of recursive subdivision. Similarly, the work of Wakamatsu et al. [START_REF] Wakamatsu | Static modeling of linear object deformation based on differential geometry[END_REF] defines another static model of elastic rods which also relies on optimization techniques. In all the cited approaches, the use of numerical methods make them too computationally expensive.

A different category of approaches, motivated by applications in computer graphics, relies on simplified deformation models that do not take into account mechanical properties. For instance, in [START_REF] Burc ¸han Bayazit | Probabilistic roadmap motion planning for deformable objects[END_REF] Bayazit et al. pre-compute some reduced models of the deformable object and use them to build a weighted roadmap through Probalistic Roadmap (PRM) methods. In [START_REF] Gayle | Path planning for deformable robots in complex environments[END_REF], Gayle et al. use the Constraint Based Motion Planning framework to simulate a deformable robot using a mass-spring model along a rough estimate of the solution path. Kabul et al. [START_REF] Kabul | Cable route planning in complex environments using constrained sampling[END_REF] extended this work to DLOs using a kinematic chain as deformable model. In these cases, the quality of the solution mostly depend on the estimate and might be unable to solve some cases.

Another interesting direction has been investigated by Mahoney et al. [START_REF] Mahoney | Deformable robot motion planning in a reduced-dimension configuration space[END_REF] by the use of a two-step process. A learning phase first collects high-dimensional samples (e.g. using simulation) and computes a new basis for the deformation set using linear dimensionality reduction. Then, this reduced deformation space is used as parameterization to perform motion planning. Its main drawback lies in the limitation to linear reduction and its nonlinear extension requires additional information on the model.

In conclusion, to overcome the curse of dimensionality, we need a model that respects the physical properties of deformations, but also minimizes the number of model parameters. This work tries to tackle all the drawbacks of the mentioned approaches. We build on top of recent work from Bretl and McCarthy [START_REF] Bretl | Quasi-static manipulation of a kirchhoff elastic rod based on a geometric analysis of equilibrium configurations[END_REF], which directly provides a parameterization of the finite-dimensional manifold of equilibrium configurations for a inextensible Kirchhoff elastic rod. This enables faster computation of deformation states and, coupled with neighborhood information, we show that it can be efficiently used for motion planning.

II. APPROACH

The work from [START_REF] Bretl | Equilibrium configurations of a kirchhoff elastic rod under quasi-static manipulation[END_REF] offers a single global chart to describe the manifold of equilibrium configurations of an elastic rod. As this chart is a low dimensional Euclidean space, it is especially well suited for sampling-based methods. However, computing the parameterization of this manifold which gives us the Direct Geometric Model (DGM) is computationally expensive with a realistic number of nodes to modelize the rod (between 100 to 300 nodes). In this case, computing time of DGM is the same order than collision checking and this constrain imply a new bottleneck for motion planning.

Our approach approximates the neighborhood of this parameterization and enables fast collision checking in this neighborhood, which is applied as an efficient approach to motion planning for free-flying quasi-static elastic rods. This approximation is based on the computation of the Jacobian at each node of the cable. Using approximate rod geometry is 10 times faster than the full DGM computation. Then, it can be efficiently used to validate a local path.

An important point is that our approximation is close to the exact solution. We show that the approximation error induced by Jacobian relative to the distance in the vicinity of the configuration space is small This low error allows us to use the approximation of the model in the exploration of the configuration space efficiently. We present our motion planning algorithm based on this Jacobian approximation.

III. MOTION PLANNING FOR QUASI-STATIC ELASTIC RODS

A. Problem Statement

Let R be a robot having n independent degrees of freedom. The configuration space C of R is an n-manifold consisting in the set of all configurations of R. Then the classic motion planning problem can be stated as finding a continuous path τ : [0, 1] → C f ree , where C f ree ⊆ C is the set of valid configurations. In general, a configuration is valid if it is collision-free in the workspace.

Consider now that the robot R is deformable, and more specifically it is an elastic Kirchhoff rod in quasi-staticity. [START_REF] Bretl | Quasi-static manipulation of a kirchhoff elastic rod based on a geometric analysis of equilibrium configurations[END_REF] showed that the set of deformations D is an m-manifold, then it is possible to apply the same formulation.

Sampling-based methods relies on a parameterization of the manifold. Unfortunatly, finding a global parameterization of a manifold is not always possible without going to higher dimensions (e.g. the Lie Group SO(3) ). As D is a manifold, it can be described by an atlas {(φ α , U α )}, that is a collection of local charts (φ α , U α ) where U α ⊆ D.

Ideally, the model of deformation should provide a single chart (φ, D), which would induce a global parameterization φ -1 of D, with φ -1 : U → D and U ⊂ R m . This parameterization would offer an elegant way to sample on the manifold and the Euclidean topology of the coordinates set gives a straightforward metric.

In this direction, the following section will introduce a global parameterization for Kirchhoff quasi-static elastic rods.

B. Kirchhoff Quasi-Static Elastic Rods

This section briefly presents the results of Bretl and McCarthy [START_REF] Bretl | Quasi-static manipulation of a kirchhoff elastic rod based on a geometric analysis of equilibrium configurations[END_REF]. We encourage the interested reader to refer to the original paper for details. Consider a Kirchhoff elastic rod with fixed end positions. At equilibrium configurations, the rod locally minimizes its total elastic energy and its shape can be expressed as a local solution to a geometric optimal control problem. From Theorem 6 in [START_REF] Bretl | Quasi-static manipulation of a kirchhoff elastic rod based on a geometric analysis of equilibrium configurations[END_REF], the set of solutions to this optimal control problem is a 6-manifold that can be parameterized by an atlas having a single chart, with coordinates given by the open subset A ⊂ R 6 . Let t ∈ [0, 1] be the parameterization along the rod and g : [0, 1] × A → SE(3) be the mapping that describes the rod spatial position at its parameter t for a given coordinate a ∈ A in previously described chart. Note these spatial positions are relative to the rod base g(0, a). For a fixed, computing g(t, a) for t ∈ [0, 1] requires the integration of several nonlinear differential systems (see Fig. 6 in [START_REF] Bretl | Quasi-static manipulation of a kirchhoff elastic rod based on a geometric analysis of equilibrium configurations[END_REF]). More specifically, sufficient conditions for static equilibrium of the rod have to be checked. From Theorem 4 in [START_REF] Bretl | Quasi-static manipulation of a kirchhoff elastic rod based on a geometric analysis of equilibrium configurations[END_REF], a rod configuration is in static equilibrium if and only if det(J(t, a)) = 0 for all t ∈ (0, 1] where the Jacobian matrix J ∈ R 6×6 is defined by:

J(t, a) = T g(t,a) L g(t,a) -1 ∂g(t, a) ∂a 1 • • • ∂g(t, a) ∂a 6 (1)
At a given rod position t, this Jacobian matrix describes variations of the geometry ∂g(t, a) with respect to variations ∂a in the chart expressed at the identity element of SE(3).

The set of coordinates a for which the rod is in static equilibrium will be noted A stable with A stable ⊂ A. These results enable us to use the coordinates A stable to describe an equilibrium configuration of the rod.

C. Planning for Free-Flying Elastic Rods

In sampling-based methods for motion planning, a key predicate is determining if a configuration lies in C f ree , i.e. checking collision in C. In most cases, this is performed by going back to the workspace W and carrying out two steps. First, we compute the DGM, which consists in a mapping from a robot configuration in C to its corresponding geometry in W. Then, we check for collisions in W between robot geometry and obstacles. If the geometry is collision free, then q lies in C f ree .

Using the configuration space formulation, adding a 6 degree of freedoms free-flyer joint to the deformable robot changes the configuration space to C = A × SE(3) which is a 12-manifold. A configuration in this space will be noted q T = a p with a ∈ A and p ∈ SE(3). Note that we consider here the free-flyer joint is attached to the base of the cable. Then the corresponding rod geometry is described by w : [0, 1] × C → W with W ⊂ SE(3). Let C col ⊆ C be the set of all configurations of the rod that lead to collisions with the environment in the workspace. Obviously, this set depends on both its deformation state a and its position p.

Let now C self = {(a, p) ∈ C | a ∈ A self }, with A self ⊆ A,
be the set of all deformation states which results in a selfcolliding geometry. Then, the set of invalid configurations is given by:

C inv = C col ∪ C self ∪ C unstable (2)
with

C unstable = {(a, p) ∈ C | a / ∈ A stable }.
Finally, we have defined the set of valid configurations for a free-flying quasi-static elastic rod, i.e. C f ree = C \ C inv .

D. Collision Checking With Fast DGM Approximation

It is now well known that the collision-checking step is a bottleneck for motion planning problems. This is especially true for virtual prototyping applications where the CAD models have generally more than 100'000 polygons. In common cases, computing the DGM is extremely fast (e.g. a kinematical chain). But in the deformable object context, computation of the DGM is considerably more costly. Using the rod model described in III-B, the DGM requires the numerical integration along the parameter t of several nonlinear differential systems. Even if these results are much more efficient than a numerical approach based on minimal energy optimization, this step becomes at least as costly as collision checking in the workspace for a sufficiently reasonable number of rod nodes. Note that the number of nodes is given by the integration resolution, and consequently the computational time to compute the DGM is in linear complexity with the number of rod nodes.

1) Fast Neighborhood Approximation: To ensure a rod configuration q is in static equilibrium, we have to compute the Jacobian J(a) of the mapping g. Without additional computational time, we can take advantage of this computation which describes the behavior of the rod geometry in the neighborhood of a in the deformation space A.

As illustrated in Fig. 1, Jacobian matrix enables us to approximate rod geometry up to 10 times faster than the full DGM computation.

2) Local path validation: In addition to check if a configuration lies in C f ree , sampling-based methods also requires a more general predicate. The local path validation consists in checking if a local path (i.e. the shortest path between two configurations with respect to the metric associated to C), lies in C f ree . This is usually done by sampling along this local path with a given resolution ∆q and checking collision in C for each sample. In our case, the geodesic between two configurations q T 1 = a 1 p 1 and q T 2 = a 2 p 2 of Fig. 2. Correspondance between approximation of the rod geometry w(q) in the workspace W (right) and approximation of the environment colliding space Ccol in the configuration space C (left). In this case, the approximation of the geometry induces a false invalid configuration detection. 3) is given by q(λ) T = a(λ) p(λ) with

C = A × SE(
a(λ) = λa 1 + (1 -λ)a 2 (3)
as A are coordinates on a chart and p(λ) is the geodesic on SE(3) as defined in [START_REF] Park | Distance metrics on the rigid-body motions with applications to mechanism design[END_REF]. The approximated geometry w(q) given by (??) can be checked for self-collisions and collisions with the environment, i.e. determining if q lies in Ccol ∪ Cself , where Ccol (resp. Cself ) are local approximations of C col (resp. C self ) due to the use of approximated geometry w(q) (see Fig. 2). Note that the locally approximated invalid configurations set Ccol ∪ Cself differs from C inv defined in (2) by not taking into account unstable rod configurations C unstable in addition to the error induced by the approximation. However, as it will be presented next, this difference is still much lower than C col ∪ C self and consequently, checking the approximated geometry w(q) is a good guess of determining if q lies in C inv .

As sampled configuration along a local path are typically close, this approximation can be efficiently used to approximately validate a local path. There are different approaches to this end, and one is illustrated in Fig. 3. The basic idea states as follow: while sampling along the local path with a resolution ∆q, we select the closest configuration where the DGM has been computed and we use its Jacobian to Fig. 3. Local path validation using DGM Jacobian approximation. The approximation of the geometry w(q i ) is obtained using the geometry w(q j ) of the closest configuration q j having exact geometry. Crosses (resp. dots) represent configurations where the DGM will be approximated (resp. fully computed).

approximate the geometry of the current sample if the distance between the two configurations is less than a threshold ∆ C ex . Otherwise, if there is no configuration with an evaluated DGM in the neighborhood bounded by ∆ C ex , a new full computation of the DGM will be performed.

This validation scheme raises two sub-problems that can be handled efficiently. First, we must choose a good value for the sampling distance ∆q. A too coarse resolution would lead to missed collisions and, in the opposite, considerable time would be wasted in unnecessary collision-checking for a too fine value. Although variations ∂w(t, q) with respect to ∂p can be computed globally, variations ∂w(t, q) with respect to ∂a require more attention as the mapping g(t, a) is only a local diffeomorphism (see proof of Theorem 8 in [START_REF] Bretl | Quasi-static manipulation of a kirchhoff elastic rod based on a geometric analysis of equilibrium configurations[END_REF]). To handle this problem, we chose to use ∆q as an initial guess of the sampling resolution and we ensure the maximum node to node distance as defined in ( 4) is less than a given maximum penetration in the workspace.

The second problem is to characterize the neighborhood where we assume the validity of our approximation. This is detailed in the following paragraph.

3) Approximation error: Our work has been highly motivated by the low distance error between an approximated rod geometry using the Jacobian and the exact geometry given by the DGM (see Fig. 4). The distance function considered here is the maximum node to node distance given by:

ρ(w 1 , w 2 ) = max t∈[0,1] d T (w 1 (t), w 2 (t)) (4)
where d T : SE(3) × SE(3) → R is the Euclidean distance between the translation part of the nodes. Thanks to this results, we have a good estimation about how far from a configuration our approximation is assumed to be valid. For a given distance error, the corresponding norm of the variation in the deformation space A denoted ∆ A ex is then used to compute its corresponding maximum variation ∆ C ex in the full configuration space C for a given local path. This enables us to take advantage of the decoupling between the two subspaces A and SE(3). This distance ∆ C ex allows us to control the size of the neighborhood we want to approximate. As presented in the following algorithm, choosing a bad value for ∆ C ex would not be critical for the validity of the result path and would only penalize performances. 4) Fast Forward Geometry RRT: FFG-RRT: In this section, we present an extension of the classical RRT algorithm which encapsulates our fast neighborhood approximation. We will call it Fast Forward Geometry RRT (FFG-RRT). The global structure of the FFG-RRT algorithm detailed in Alg. 1 is very similar to a Lazy-RRT. In our context of motion planning for quasi-static Kirchhoff elastic rods, the APPROX EXTEND function could be the local path validation as described in III-D.2.

The second part of the algorithm is given by the CHECK SOLUTION function. This function tries to find a solution path in the current tree where most of the edges are only approximately valid and check for exact validity of this path. Invalidated edges are then removed from the tree.

Algorithm 1 FFG-RRT(q start , q goal ) 1: T .init(q start ) 2: for i ← 1 to k do 3:

q rand ← RANDOM SAMPLE() 4:
q near ← NEAREST(T , q rand ) 5:

q new ← APPROX EXTEND(q near , q rand ) 6:

if q new = q near then 7:

T .addVertex(q new ) 8:

T .addEdge(q near , q new ) 9: end if 10:

(solved, τ ) ← CHECK SOLUTION(q start , q goal , T ) 11:

if solved then return τ 12:

end if 13: end for 14: return failure Note the efficiency of this algorithm relies on the number of tree reconstructions in the function CHECK SOLUTION. As we seen that our approximation is close to the exact DGM, this implies a low number of calls to this function and this will be verified in the experimental results. Fig. 6. Two industrial cases of disassembling for a free-flying cable (top) and a fixed base cable (bottom). The cable has to get out from a highly constrained start configuration (in green) to a low constrained configuration (in red). Closer views on the start configuration are represented on the right for each scenario.

• Fixed engine. The model is the same as in the previous scenario but here the cable has a fixed base. Note that on all of these scenarios, the rod has to deform its shape to go from start to goal configurations. We ran the benchmark 30 times for each case. The number of rod nodes has been chosen to give sufficient realism of the physical deformation and depends on the level of deformation and the rod length. Benchmark results are shown in Table I. Times are average planning time in seconds and success rate is in percentage. We see our approach solves about two to three times faster than its respective classical implementation, with a similar success rate. Also, as toy scenarios typically show two connected components with a narrow passage, bi-directional planners such RRT-Connect are highly more efficient in these cases. Thanks to the quality of the approximation, the function CHECK SOLUTION is rarely called (0 to 2 calls) which ensures us a significant performance improvement. Theses cases where the approximated solution path is invalidated are mainly due to the non respect of sufficient conditions for static equilibrium as described in III-B (i.e. q ∈ C unstable ).

C. Results

V. CONCLUSION AND PERSPECTIVES

In this paper, we presented an new approach to the motion planning problem for a Kirchhoff elastic rod. Using a realistic physical model for the DLO, we are able to solve complex industrial scenarios in a reasonable time. We also introduced a new local planning scheme by taking advantage of the linearization of the Direct Geometrical Model which enables to solve the problem two to three times faster. This local method can be embedded in many different motion planning algorithms and extends to any model where the DGM cost becomes critical.

In the future, we would like to extend the approach by coupling with dynamic simulation to ensure mechanical constraints of the rod between two quasi-static rod configurations. Another interesting problem in manufacturing would be to take into account gripers of the rod during the planning process.

Fig. 1 .

 1 Fig. 1. Log scale computation time of the DGM (red) and its first order approximation (green).

Fig. 4 .

 4 Fig. 4. Measurements of the distance error between an approximated and an exact geometry depending on the norm of the variation δa in A.

ACKNOWLEDGMENT

This work was supported by the French National Research Agency under the project Flecto (ANR-Digital Models). Industrial models are courtesy of Siemens-KineoCAM.

The two toy scenarios used for benchmarks: crack (left) and backward (right). Start (resp. goal) configurations of the rod are represented in green (resp. red). Note the head of cable has been marked with a magenta circle on the Benchmark scenario to emphasize solution path could not be obtained with a point-like robot.

IV. EXPERIMENTAL RESULTS

In this section we will present and analyze results on solving the motion planning problem for free flying Kirchhoff elastic rods using the model presented in section III-B and our approach for fast approximation of the Direct Geometric Model detailed in section III-D.

A. Implementation Detail

We chose to implement our approach on top of the classical RRT and RRT-Connect [START_REF] Kuffner | Rrt-connect: An efficient approach to single-query path planning[END_REF] planners for the benchmarks using the C++ motion planning library OMPL [START_REF] Ioan | The Open Motion Planning Library[END_REF]. The numerical integration of differential systems required to obtain the geometrical state of the rod was implemented using fourth order Runge-Kutta methods. For collision checking, the geometry of the rod was approximated to a hierarchical chain of capsules (i.e. Line Swept Spheres). The Flexible Collision Library (FCL) [START_REF] Pan | Fcl: A general purpose library for collision and proximity queries[END_REF] was used to perform collision checking computations. All the benchmarks were run on a PC with 8GB of main memory and using one core of an Intel Core i7-2720QM processor running at 2.2Ghz.

B. Benchmark Scenarios

We selected four distinguished scenarios in order to test the effectiveness of our approach. Two firsts consist in toy scenarios (see Fig. 5) where each shows a specific difficulty. Last two scenarios are industrial cases with a disassembly study (see Fig. 6). On all of this cases, bounds have been set on the configuration space but, for clarity, environment bounding boxes are not shown in the illustrations.

• Crack. Lightweight model where the rod must pass through a crack shaped narrow passage. The length (resp. diameter) of the rod is half the length (resp. width) of the crack free space. • Backward. Lightweight model where a point-like robot path would give an invalid guess or approximate solution. Note that due to the bounding box, the DLO has to go trough the corridor between obstacles. • Free-Flying engine. Industrial model with 132,000 polygons and where the rod models a cable in a typical disassembly study case.


Bartłomiej Błaszczyszyn
email: bartek.blaszczyszyn@ens.fr

Miodrag Jovanovic
email: miodrag.jovanovic@orange.com

Mohamed Kadhem Karray
email: mohamed.karray@orange.com

How user throughput depends on the traffic demand in large cellular networks

Keywords: user-througput, traffic demand, cell-load, Little's law, cellular network, typical cell, point process, ergodicity, Palm theory, measurements

We assume a space-time Poisson process of call arrivals on the infinite plane, independently marked by data volumes and served by a cellular network modeled by an infinite ergodic point process of base stations. Each point of this point process represents the location of a base station that applies a processor sharing policy to serve users arriving in its vicinity, modeled by the Voronoi cell, possibly perturbed by some random signal propagation effects. User service rates depend on their signal-to-interference-and-noise ratios with respect to the serving station.

Little's law allows to express the mean user throughput in any region of this network model as the ratio of the mean traffic demand to the steady-state mean number of users in this region. Using ergodic arguments and the Palm theoretic formalism, we define a global mean user throughput in the cellular network and prove that it is equal to the ratio of mean traffic demand to the mean number of users in the steady state of the "typical cell" of the network. Here, both means account for double averaging: over time and network geometry, and can be related to the per-surface traffic demand, base-station density and the spatial distribution of the signal-to-interference-and-noise ratio. This latter accounts for network irregularities, shadowing and cell dependence via some cell-load equations.

Inspired by the analysis of the typical cell, we propose also a simpler, approximate, but fully analytic approach, called the mean cell approach. The key quantity explicitly calculated in this approach is the cell load. In analogy to the load factor of the (classical) M/G/1 processor sharing queue, it characterizes the stability condition, mean number of users and the mean user throughput. We validate our approach comparing analytical and simulation results for Poisson network model to real-network measurements.

I. INTRODUCTION

Mean user throughput is a key quality-of-service metric in cellular data networks. It describes the average "speed" of data transfer during a typical data connection. It is usually defined as the ratio of the average number of bits sent (or received) per data request to the average duration of the data transfer. Since coexisting connections in a given network cell share some given cell transmission capacity, mean user throughput depends inherently on the requested data traffic. It also depends on the network architecture (positioning of the base stations) and in fact may significantly vary across different network cells. Moreover, extra-cell interference makes performance of different cells interdependent. Predicting the mean user throughput in function of the mean traffic demand locally (for each cell) and globally in the network (which involves appropriate spatial averaging in conjunction with the temporal one, already present in the classical definition of the throughput) is a key engineering task in cellular communications. In this paper we propose an analytic approach to the evaluation of the mean user throughput in large irregular cellular networks, validated by real-network measurements performed in operational networks.

Little's law allows to calculate the mean user throughput as the ratio of the mean traffic demand (number of bits requested per unit of time) to the mean number of users in the steady state of the network. This argument can be used to express mean user throughput locally in any region of the network. Statistics usually collected in operational networks allow to estimate the mean traffic demand and the mean number of users for each cell and hour of the day. Even if they carry important information about the local network performance, they exhibit important variability over time (24 hours) and network cells; cf. Figure 1. This can be explained by the fact that mean user throughput in a particular cell does not depend only on the traffic in this cell, but also on the neighbouring cells. Moreover, the geometry of different cells in a real network may significantly differ. For these two reasons, the family of local (established for each cell) throughputversus-traffic laws usually exhibits a lot of variability both in real data and in network simulations, and hence does not explain well the macroscopic (network-level) relation between the mean traffic demand and mean user throughput. Finding such a macroscopic relation is an important task for network dimensioning. It is clear that an appropriate spatial averaging is necessary to discover such a macroscopic law.

Spatial averages of point patterns, modeling in our case the geographic locations of base stations, can be studied using the formalism of Palm distributions naturally related to the ergodic results for point processes. Within this setting one considers a typical base station with its typical cell (zone of service) whose probabilistic characteristics correspond to the aforementioned spatial averages of the characteristics for all base stations in the network. Adopting this formalism, we define the mean user throughput in the infinite ergodic network as the limit of the ratio of the mean volume of the data request to the mean service duration in a large, increasing to the whole plane, network window. As the main result, we prove that such defined (macroscopic) throughput characteristic is equal to the ratio of the mean traffic demand to the mean number of users in the typical cell of the network. Both these means account for double averaging: over time and network geometry.

A key element of the analysis of the cellular network is the spatial distribution of the signal-to-interference-andnoise ratio (SINR). We show how this distribution enters into the macroscopic characterization of the throughput. When considering SINR we are able to account for the fact that the base stations that are idling, i.e., have no users to serve, do not contribute to the interference. This makes the performance of different cells interdependent and we take it into account via a system of cell-load equations.

Finally, we show how to amend the model letting it account for the shadowing in the path loss. The latter is known to impact the geometry of the network, in the sense that the serving base station is not necessarily the closest one. It also alters the distribution of the SINR.

A. Related work

The evaluation of user QoS metrics in cellular networks is a hard problem, but crucial for network operators and equipment manufacturers. It motivates a lot of engineering and research studies. The complexity of this problem made many actors develop complex and time consuming simulation tools such as those developed by the industrial contributors to 3GPP (3rd Generation Partnership Project) [START_REF] Gpp | TR 36.814-V900 Further advancements for E-UTRA -Physical Layer Aspects[END_REF]. There are many other simulation tools such as TelematicsLab LTE-Sim [START_REF] Piro | Simulating LTE cellular systems: an open source framework[END_REF], University of Vien LTE simulator [START_REF] Mehlführer | The Vienna LTE simulators -Enabling reproducibility in wireless communications research[END_REF][START_REF] Simko | Optimal pilot symbol power allocation under time-variant channels[END_REF] and LENA tool [START_REF] Baldo | An open source product-oriented LTE network simulator based on ns-3[END_REF][START_REF] Baldo | A new model for the simulation of the LTE-EPC data plane[END_REF] of CTTC, which are not necessarily compliant with 3GPP.

A possible analytical approach to this problem is based on the information theoretic characterization of the individual link performance; cf e.g. [START_REF] Goldsmith | Variable-rate variablepower MQAM for fading channels[END_REF][START_REF] Mogensen | LTE Capacity Compared to the Shannon Bound[END_REF], in conjunction with a queueing theoretic modeling and analysis of the user traffic; cf. e.g. [START_REF] Borst | User-level performance of channel-aware scheduling algorithms in wireless data networks[END_REF][START_REF] Bonald | Wireless downlink data channels: user performance and cell dimensioning[END_REF][START_REF] Hegde | Capacity of multiservice WCDMA Networks with variable GoS[END_REF][START_REF] Bonald | Flow-level performance and capacity of wireless networks with user mobility[END_REF][START_REF] Rong | Performance evaluation of cellular networks offering TV services[END_REF][START_REF] Karray | A queueing theoretic approach to the dimensioning of wireless cellular networks serving variable bit-rate calls[END_REF]. All these works consider some particular aspects of the network and none of them considers a large, irregular multicell network. Such a scenario is studied in our approach by using stochastic-geometric tools combined with the two aforementioned theories. As a result, we propose a global, macroscopic approach to the evaluation of the user QoS metrics in cellular networks, which we compare and validate with respect to real network measurements.

Stochastic geometry has already been shown to give analytically tractable models of cellular networks, see e.g. [START_REF] Andrews | A tractable approach to coverage and rate in cellular networks[END_REF]. However, to the best of our knowledge, the prior works in this context usually do not consider any dynamic user traffic.

B. Paper organization

We describe our general cellular network model, comprising the geometry of base stations, path loss model, space-time traffic demand process, and the service policy in Section II. In Section III we develop two approaches, called respectively typical and mean cell approach, allowing to study the dependence of the mean user throughput in the network and other macroscopic network characteristics on the traffic demand and other model parameters. In Section IV we apply these approaches studying some particular network model and compare the obtained numerical results to real field measurements.

II. MODEL DESCRIPTION AND LOCAL

CHARACTERISTICS

We shall now describe our model.

A. Network model

We consider locations {X 1 , X 2 , . . .} of base stations (BS) on the plane R2 as a realisation of a point process, which we denote by Φ. 1 We assume that Φ is stationary and ergodic with positive, finite intensity (mean number of BS per unit of surface) λ BS . 2 In order to simplify the presentation, we shall make first the following two assumptions, which will be relaxed in Sections III-E and III-D, respectively. 1) There is no shadowing. The (time-averaged over fading)

propagation loss depends only on the distance r between the emitter and the receiver through a path-loss function l(r), which we assume increasing. 2) Full interference. Each base station is always transmitting at some fixed power P , common for all stations.

We will also assume throughout the whole paper that each user is served by the BS which he or she receives with the strongest signal power. The consequence of the assumption 1 above is that each BS X ∈ Φ serves users in a geographic zone V (X) = {y ∈ R 2 : |y -X| ≤ min Z∈Φ |y -Z|} which is called Voronoi cell of X in Φ. Both the above assumptions 1 and 2 allow to represent the signal to interference and noise ratio (SINR) at location y ∈ V (X), X ∈ Φ as SINR (y, Φ) := P/l (|y -X|)

N + P Z∈Φ\{X} 1/l (|y -Z|) , ( 1
)
where N is the noise power. Note that this represents the SINR in the down-link (BS to user) channel. Note also, that we assume that each interfering base station always transmits (even if, for example, it has no user to serve). This model, which we call full interference model, will be improved in Section III-D. We assume that the peak bit-rate at location y, defined as the bit-rate of a user located at y when served alone by its BS, is some function R (SINR), of the SINR. 3We consider variable bit-rate (VBR) traffic; i.e., users arrive to the network and require to transmit some volume of data at a bit-rate decided by the network. Specifically, we assume that user channels are intra-cell orthogonal and inter-cell independent: if BS X serves n users located at y 1 , y 2 , . . . , y n ∈ V (X) then the bit-rate of the user located at y j equals to 1/n th of its peak bit-rate 1 n R (SINR (y j , Φ)), j ∈ {1, 2, . . . , n}. 4The pattern of BS Φ does not evolve in time. We describe now the space-time process of user arrivals (and departures).

B. Traffic demand

Users arrive uniformly on the plane and require to transmit a random (arbitrarily distributed) volume of data of mean 1/µ bits. The duration between the arrivals of two successive users in each geographic zone S of surface |S| is an exponential random variable of parameter λ × |S| . This means that on average there are λ arrivals per surface unit. The arrival locations, inter-arrival durations as well as the data volumes are assumed independent across users. The above description corresponds to a space-time Poisson process of arrivals, independently marked by data volumes.

We assume that the users don't move considerably during their calls. 5 Each user stays in the system for the time necessary to download his data. This takes a random (service) time because the bit-rate with which he is served depends on the configuration of other users served by the same base station. Users depart from the system immediately after having downloaded their data.

The traffic demand per surface unit is then equal to ρ = λ µ , which will be expressed in bit/s/km 2 . The traffic demand in a given cell equals

ρ (X) = ρ |V (X)| , X ∈ Φ. (2)

C. Local quality of service characteristics

For a fixed configuration of BS Φ, the service of users arriving to the cell V (X) of a given BS X ∈ Φ can be modeled by an appropriate (spatial) multi-class processor sharing queue, with classes corresponding to different peak bitrates characterized by user locations y ∈ V (X). Note also that a consequence of our model assumptions (in particular the full interference assumption 2, inter-cell channel independence and space-time Poisson arrivals) the service processes of different queues are independent.

We consider now the steady-state of users served in each cell V (X). 6 The following expressions follow from the queueingtheoretic analysis of the processor sharing systems of each BS X ∈ Φ, cf [START_REF] Bonald | Wireless downlink data channels: user performance and cell dimensioning[END_REF][START_REF] Karray | A queueing theoretic approach to the dimensioning of wireless cellular networks serving variable bit-rate calls[END_REF] for the details.

• The service process of BS X ∈ Φ is stable if and only if its traffic demand does not exceed the critical value that is the harmonic mean of the peak bit-rate over the cell:

ρ c (X) := |V (X)| V (X) 1/R (SINR (y, Φ)) dy . (3)
Note that ρ c (X) depends on X and on Φ. The same observation is valid for the subsequent cell characteristics.

• The mean user throughput in the given cell, defined as the ratio of the mean volume of the data request 1/µ to the average service time of users in this cell, can be expressed as follows

r (X) = max(ρ c (X) -ρ (X) , 0) . (4)
• The mean number of users in the steady state of the given cell equals to

N (X) = ρ (X) r (X) . (5) Note that N (X) = ∞ if ρ(X) ≥ ρ c (X).
• The probability that the given BS is not idling in the steady state (has at least one user to serve) equals

p (X) = min (θ (X) , 1) , (6)
where θ(X), which we call cell load, is defined as

θ (X) := ρ (X) ρ c (X) . (7)
Note that the cell is stable if and only if θ(X) < 1 and

θ(X) = ρ V (X) 1/R (SINR (y, Φ)) dy . (8)
Moreover,

N (X) = θ(X) 1 -θ(X) , (9) r
(X) = ρ(X)(1/θ(X) -1) (10)
provided θ(X) < 1.

The above expressions allow to express all other characteristics in terms of the traffic demand per cell ρ(X) and the cell load θ(X).

Remark 1: All the above characteristics are local network characteristics in the sense that they characterize the service at each BS X and vary over X ∈ Φ. Real data analysis and simulations for Poisson network models exhibit a lot of variability among these characteristics. In particular, plotting the mean user throughput r(X) as function of the mean traffic demand ρ(X) for different X ∈ Φ does not reveal any apparent systematic relation between these two local characteristics; cf. Figure 1.

III. GLOBAL NETWORK CHARACTERISTICS

In this section we propose some global characteristics of the network allowing to characterize its macroscopic performance. We are particularly interested in finding such a relation between the (per surface) traffic demand ρ and the (global) mean user throughput in the network, with this latter characteristic yet to be properly defined.

A. Typical cell of the network

A first, natural idea in this regard is to consider spatial averages of the local characteristics in an increasing network window A, say a ball centered at the origin and the radius increasing to infinity. Assuming ergodicity of the point process Φ of the BS, these averages can be expressed and calculated as Palm-expectations of the respective characteristics of the so called "typical cell" V (0). For example

lim |A|→∞ 1/Φ(A) X∈A ρ(X) = E 0 [ρ(0)] = ρE 0 [|V (0)|] . (11)
The typical cell V (0) is the cell of the BS located at the origin X = 0 and being part of the network Φ distributed according to the Palm distribution P 0 associated to the original stationary distribution P of Φ. In the case of Poisson process, the relation between the Palm and stationary distribution is particularly simple and (according to Slivnyak's theorem) consists just in adding the point X = 0 to the stationary pattern Φ.

The convergence analogue to [START_REF] Hegde | Capacity of multiservice WCDMA Networks with variable GoS[END_REF] holds for each of the pre-

viously considered local characteristics E 0 [ρ c (0)], E 0 [r(0)], E 0 [N (0)], E 0 [p(0)] and E 0 [θ(0)].
The convergence is P almost sure and follows from the ergodic theorem for point processes (see [START_REF] Baccelli | Stochastic Geometry and Wireless Networks, Volume I -Theory, ser. Foundations and Trends in Networking[END_REF]Theorem 4.2.1], [START_REF] Daley | An introduction to the theory of point processes[END_REF]Theorem 13.4.III]). However, as we will explain in what follows, not all of these mean-typical cell characteristics have natural interpretations as macroscopic network characteristics.

First, note that the existence of some (even arbitrarily small) fraction of BS X which are not stable (with ρ(X)

≥ ρ c (X), hence N (X) = ∞) makes E 0 [N (0)] = ∞.
Remark 2: For a well dimensioned network one does not expect unstable cells. For a perfect hexagonal network model Φ all cells are stable or unstable depending on the value of the per-surface traffic demand ρ. An artifact of an infinite, homogeneous, Poisson model Φ is that for arbitrarily small ρ there exists a non-zero fraction of BS X ∈ Φ, which are nonstable. This fraction is very small for reasonable ρ, allowing to use Poisson to study QoS metrics which, unlike E 0 [N (0)], are not "sensitive" to this artifact.

We will also show in the next section that it is not natural to interpret E 0 [r(0)] (which is not sensitive to the existence of a small fraction of unstable cells) as the mean user throughput in the network; see Remark 7. Before we give an alternative definition of this latter QoS, let us state the following result, which will be useful in what follows.

Proposition 3: We have

E 0 [ρ(0)] = ρ λ BS , (12)
E 0 [θ(0)] = ρ λ BS E[1/R (SINR (0, Φ))] . (13)
Proof: The first equation is quite intuitive: the average cell surface is equal to the inverse of the average number of BS per unit of surface. Formally, both equations follow from the inverse formula of Palm calculus [17, Theorem 4.2.1]. In particular, for (13) one uses representation (8) in conjunction with the inverse formula.

Remark 4: Note that the expectation in the right-hand-side of ( 13) is taken with respect to the stationary distribution of the BS process Φ. It corresponds to the spatial average of the inverse of the peak bit-rate calculated throughout the network. The (only) random variable in this expression is the SINR experienced by the typical user. This distribution is usually known in operational networks (estimated from user measurements). It can be also well approximated using Poisson network model for which its distribution function admits an explicit expression; cf [START_REF] Błaszczyszyn | Using Poisson processes to model lattice cellular networks[END_REF][START_REF] Keeler | SINRbased k-coverage probability in cellular networks with arbitrary shadowing[END_REF].

B. Mean user throughput in the network

Faithful to the usual definition of the mean user throughput as the ratio of the mean volume of the data request to the mean service duration (which we retained at the local, cell level) we aim to define now the mean user throughput in the (whole) network as the ratio of these two quantities taken for increasing network window A. However in order to "filter out" the impact of cells which are not stable and avoid undesired degeneration of this characteristic (e.g. for Poisson process; cf. Remark 2) let us consider the union of all stable cells

S := X∈Φ:ρ(X)<ρc(X) V (X) .
Note that the stationarity of Φ implies the same for the random set S. We denote by π S = E[1(0 ∈ S)] the volume fraction of S and call it the stable fraction of the network. It is equal to the average fraction of the plane covered by the stable cells; cf. [17, Definition 3.4 and the subsequent Remark]. Denote also

N 0 := E 0 [N (0)1(N (0) < ∞)] . (14)
We are ready now to define the mean user throughput in the network r 0 as the ratio of the average number of bits per data request to the average duration of the data transfer in the stable part of the network

r 0 := lim |A|→∞ 1/µ (temporal-)mean service time in A ∩ S . (15)
Here is the key result of the typical cell approach. Its proof is given at the end of this section. Theorem 5: For ergodic network Φ we have

r 0 = ρ π S λ BS N 0 . ( 16
)
Remark 6: Equation ( 16) provides a macroscopic relation between the traffic demand and the mean user throughput in the network, which we are primarily looking for in this paper. It will be validated by comparison to real data measurements. Quantities N 0 and π S do not have explicit analytic expressions analogous to [START_REF] Rong | Performance evaluation of cellular networks offering TV services[END_REF]. Nevertheless they can be estimated from simulations of a given network model Φ. Note that these are static simulations of the network model. No simulation of the traffic demand process is necessary, which greatly simplifies the task. For small and moderate values of the traffic demand (observed in real networks) one obtains π S 1. Moreover, in Section III-C we will propose some more explicit approximation of N 0 .

Remark 7: Assume that there are no unstable cells in the network. This is the case e.g. for lattice (say hexagonal) network models with traffic demand ρ < ρ c (X) = ρ c , where the value of the critical traffic is the same for all cells. Then π S = 1, N 0 = E 0 [N (0)] and the relation ( 16) takes form

r 0 = ρ λ BS E 0 [N (0)] = E 0 [ρ(0)] E 0 [N (0)] . ( 17
)
Thus, in general

r 0 = E 0 [r(0)] = E 0 [ρ(0)/N (0)].
We want to emphasize that this is not merely a theoretical detail resulting from our (and common) definition of the mean throughput [START_REF] Andrews | A tractable approach to coverage and rate in cellular networks[END_REF]. The expression E 0 [r(0)] = E 0 [ρ(0)/N (0)], which in principle can be considered as another global QoS metric, is in practice difficult to estimate. Indeed, when estimating E 0 [r(0)] as the average of the ratio "traffic demand to the number of users" from real data measurements, one needs to give a special treatment to observations which correspond to cells during their idling hours (i.e., with no user, and such observations are not rare in operational networks). Neither skipping nor literal acceptance of these observations captures the right dependence of the mean user throughput on the traffic demand.

Proof of Theorem 5: By Little's law [20, eq. (3.1.14)] the temporal mean service time T W of users in any region of the network W , say the union of stable cells with BS in some region A, W = X∈A∩S V (X), is related to the mean number N W of the users served in this region W in the steady state by the equation N W = λ|W |T W . Consequently, the mean user throughput in this region W can be expressed as

1/(µT W ) = ρ|W |/N W . Using |W | N W = |W | X∈A∩S N (X) (18) = X∈A |V (X)|1(N (X) < ∞) |A| |A| X∈A∩S N (X) (19)
and again the ergodic theorem for point yprocess Φ, we obtain the that limit in ( 15) is P-almost surely equal to

ρE 0 [|V (0)|1(N (0) < ∞)]/E 0 [N (0)1(N (0) < ∞)]
. By the aforementioned inverse formula of Palm calculus we conclude

E 0 [|V (0)|1(N (0) < ∞)] = E[1(0 ∈ S)]/λ BS .

C. Mean cell

It is tempting to look for a synthetic model which would allow to relate main parameters and QoS metrics of a large irregular cellular network in a simple, yet not simplistic way. The typical cell approach described up to now offers such possibility. In this section we will go a little bit further and propose an even simpler model. It consists in considering a virtual cell, to which we will assign the parameters and QoS metrics inspired by the analysis of the typical cell. In contrast to the typical cell, our virtual cell is not random and this is why we call it the mean cell. Specifically, we define it as a (virtual) cell having the same traffic demand ρ and load θ as the typical cell. Note that these two characteristics admit explicit expressions; cf. Proposition 3 and Remark 4.

ρ := E 0 [ρ (0)] = ρ λ BS , (20)
θ := E 0 [θ (0)] = ρ λ BS E[1/R (SINR (0, Φ))] . (21)
For the remaining characteristics, we assume that they are related to the above two via the relations presented in Section II-C. Specifically, following (2) we define the surface of the mean cell by V = ρ/ρ and in analogy to [START_REF] Goldsmith | Variable-rate variablepower MQAM for fading channels[END_REF] we define the critical load of the mean cell as ρc := ρ θ .

We say that the mean cell is stable if ρ < ρc . Inspired by (4) we define the user's throughput in the mean cell by r := max (ρ c -ρ, 0)

and, as in ( 5), the mean number of users in the mean cell is defined as N := ρ r .

We observe the following immediate relations. Corollary 8: The mean cell is stable if and only if θ < 1. In this case

N = θ 1 - θ , (23)
r = ρ(1/ θ -1) , (24)
which are analogous to ( 9) and ( 10), respectively. Remark 9: The equation ( 24) provides an alternative macroscopic relation between the traffic demand and the mean user throughput in the network. It is purely analytic; no simulations are required provided one knows the distribution of the SINR of the typical user in [START_REF] Siomina | Analysis of cell load coupling for lte network planning and optimization[END_REF]. It will be validated by comparison to real data measurements. We consider it as an approximation of [START_REF] Daley | An introduction to the theory of point processes[END_REF]. It consists in assuming N N 0 /π S . This latter hypothesis will be also separately validated numerically.

Remark 10: Note that the key characteristic of the mean cell is its load θ. In analogy to the load factor of the (classical) M/G/1 processor sharing queue, it characterizes the stability condition, mean number of users and the mean user throughput.

D. Cell-load equations

We have to revoke now the full interference assumption 2 made in Section II-A. An amendment is necessary in this matter for the model to be able to predict the real network data; cf numerical examples in Section IV. Recall that the consequence of this assumption is that in the expression (1) of the SINR all the interfering BS are always transmitting at a given power P . In real networks BS transmit only when they serve at least one user. 7 Taking this fact into account in an exact way requires introducing in the denominator of (1) the indicators that a given station Z ∈ Φ at a given time is not idling. This, in consequence, would lead to the probabilistic dependence of the service process at different cells and result in a non-tractable model. In particular, we are not aware of any result regarding the stability of such a family of dependent queues. For this reason, we take into account whether Z is idling or not in a simpler way, multiplying its powers P by the probability p(Z) that it is not idle in the steady state. In other words we modify the expression of the SINR as follows SINR (y, Φ) := P/l (|y -X|)

N + P Z∈Φ\{X} p (Z) /l (|y -Z|) , (25)
for y ∈ V (X), X ∈ Φ where p(Z) are cell non-idling probabilities given by ( 6). We will see in Section IV that this model, called (load-)weighted interference model, fits better to real field measurements than the full interference model. The above modification of the model preserves the independence of the processor-sharing queues at different cells given the realization Φ of the network (thus allowing for the explicit analysis of Section II-C). However the cell loads θ(X) are no longer functions of the traffic demand and the SINR experienced in the respective cells, but are related to each other by the following equations that replace ( 8)

θ (X) = ρ V (X) 1 R P/l(|y-X|) N +P Z∈Φ\{X} min(θ(Z),1)/l(|y-Z|)
dy .

(26) We call this system of equations in the unknown cell loads {θ (X)} X∈Φ the cell-load equations.

Remark 11: [Spatial stability] The weighted interference model introduces more "spatial" dependence between the processor sharing queues of different cells, while preserving their "temporal" (conditionally, given Φ) independence. A natural question regarding the existence and uniqueness of the solution of the fixed point problem (26) arises. Note that the mapping in the right-hand-side of (26) is increasing in all θ(Z), Z ∈ Φ provided function R is increasing. Using this property it is easy to see that successive iterations of this mapping started off θ(Z) ≡ 0 on one hand side and off θ(Z) as in [START_REF] Mogensen | LTE Capacity Compared to the Shannon Bound[END_REF] (full interference model) on the other side, converge to a minimal and maximal solution of (26), respectively. An interesting theoretical question regards the uniqueness of the solution of (26), in particular for a random, say Poisson, point process Φ. Answering this question, which we call "spatial stability" of the model, is unfortunately beyond the scope of this paper. 8 The simulation study of the typical cell model, presented in Section IV (where we use Matlab to find a solution of (26) for any given finite pattern of base stations Φ) is less stable for larger values of the traffic demand ρ.

In the mean cell approach (cf Section III-C) we take into account the weighted interference model by the following (single) equation in the mean-cell load

θ θ = ρ λ BS E 1/R P/l (|X * |) N + P Z∈Φ\{X * } θ/l (|y -Z|) , (27)
where X * is the location of the BS whose cell covers the origin. We solve the above equation with θ as unknown. We will show in the numerical section that the solution of this equation gives a good estimate of the empirical average of the 8 Existence and uniqueness of the solution of a very similar problem (with finite number of stations and a discrete traffic demand) is proved in [START_REF] Siomina | Analysis of cell load coupling for lte network planning and optimization[END_REF].

loads {θ (X)} X∈Φ obtained by solving the system of cell-load equations ( 26) for the simulated model.

Remark 12: [Pilot channel] The cells which are not idle might still emit some power (e.g. in the pilot channel). This can be taken into account by replacing p(Z) = min(θ(Z), 1) in (26) by p(Z)(1 -) + , where is the fraction of the power emitted all the time. Similar modification concerns θ in the right-hand-side of ( 27).

E. Shadowing

Until now we were assuming that the propagation loss is only induced by the distance between the transmitter and the receiver. In this section we will briefly explain how the effect of shadowing can be taken into account.

Assume that the shadowing between a given station X ∈ Φ and all locations y ∈ R 2 is modeled by some random field S X (y -X). That is, we assume the propagation loss between X and y L X (y) = l(|y-X|) S X (y-X) . We assume that, given Φ, the random fields S X (•) are independent across X ∈ Φ and identically distributed. In general, we do not need to assume any particular distribution for S X (•) (neither independence nor the same distribution of S X (y) across y).

The assumption that each user is served by the BS received with the smallest path-loss results in the following modification of the geographic service zone of X, which we keep calling "cell"

V (X) = {y ∈ R 2 : L X (y) ≤ min Z∈Φ L Z (y)} . ( 28
)
For mathematical consistence we shall assume that, almost surely, the origin belongs to a unique cell (i.e., is not located on any cell boundary).

The SINR at location y can be expressed by ( 1) or (25) with l(|y -Y |) replaced by L Y (y) for Y ∈ Φ, depending whether we consider the full interference or the weighted interference model, for y ∈ V (X), with V (X) defined by (28). The same modification regards the cell-load equations ( 26) and ( 27).

All the previous results involving the typical cell remain valid for this modification of the model. In particular, the results of Proposition 3 can be extended to the model with shadowing (where the cell associated to each base station is not necessarily the Voronoi cell) provided the origin 0 belongs to a unique cell almost surely.

Note that the mean cell surface E 0 [|V (0)|] = 1/λ BS , and hence the mean traffic demand per cell E 0 [ρ(0)] = ρ/λ BS , do not depend at all on the shadowing. The values of other characteristics of the typical and the mean cell will change depending on the distribution of the random shadowing field S X (y) (both the marginal distributions and the correlation across y). An interesting remark in this regard is as follows.

In the full interference model, the mean load of the typical cell and the load of the mean cell (which are by the definition equal E 0 [θ(0)] = θ) depend only on the stationary marginal distribution of SINR(0, Φ), cf. [START_REF] Siomina | Analysis of cell load coupling for lte network planning and optimization[END_REF]. Hence, it does not depend on the (spatial) correlation of S X (y) across y. Moreover, this distribution is known in the case of the Poisson network and identically distributed marginal shadowing S X (y) ∼ S. As explained in [START_REF] Keeler | SINRbased k-coverage probability in cellular networks with arbitrary shadowing[END_REF], in this case SINR(0, Φ) in the model with shadowing has the same distribution as in the model without shadowing and the density of stations equal to

λ BS × E[S 2/β ].
In particular, a specific distribution of S and the correlation of S X (y) across y play no role. This equivalence of the two models (with and without shadowing) is more general, as explained in [START_REF] Keeler | SINRbased k-coverage probability in cellular networks with arbitrary shadowing[END_REF][START_REF] Błaszczyszyn | Equivalence and comparison of heterogeneous cellular networks[END_REF], and applies also to the mean-cell-load equation ( 27).

Remark 13: The impact of the shadowing on the mean cell model, both in full and weighted interference scenario in the above Poisson model with S X (y) ∼ S can be summarized as follows. It modifies only the cell load θ and not the traffic demand ρ. Moreover, multiplying λ BS by E[S 2/β ] and dividing ρ by the same moment one obtains an equivalent (in terms of all considered characteristics) mean cell without shadowing.

IV. NUMERICAL RESULTS

To illustrate the motivation of this work, we present first on Figure 1 non-averaged data obtained from the measurements performed in an operational network in some zone of some big city in Europe. 9 Different points in this figure correspond to the measurements of the traffic demand and the estimation of the user throughput made by different cells during different hours of the day. No apparent relation between these two quantities can be observed in this way.

In order to uderstand and predict the performance of the network for which we have presented the above data, we will now specify correspondingly our general model and study it using the proposed approach. The obtained results 9 More precisely, a dense urban network zone consisting of 382 base stations was selected in a big European city, whose locations loosely satisfy homogeneous spatial Poisson assumption. Ripley's L-function, cf [23, page 50], plotted on Figure 9, was used to verify this latter assumption. The density of base stations in this dense urban zone is about 4.62 base stations per km 2 . Later, we will consider also a urban zone of a different European city, where the spatial homogeneous Poissonianity of the base station locations can also be retained; cf. Figure 9, with roughly four times smaller density of base stations, more precisely 1.15 stations per km 2 . In both cases the network operates HSDPA system with MMSE coding.

will be compared to the appropriately averaged real field measurements.

A. Model specification

Consider the following numerical setup. Assume Poisson process of BS with intenisty λ BS = 4.62km -2 (which corresponds to an average distance between two neighbouring BS of 0.5km). We assume the path-loss function l(r) = (Kr) β , with K = 7117km -1 , and the path loss exponent β = 3.8. The propagation model comprises the log-normal shadowing with the logarithmic standard deviation 10dB; cf [START_REF] Błaszczyszyn | Using Poisson processes to model lattice cellular networks[END_REF], and the mean spatial correlation distance 0.05km.

The transmision power is P = 58dBm, with the fraction = 10% used in the pilot channel. The antenna pattern is described in [START_REF] Gpp | TR 36.814-V900 Further advancements for E-UTRA -Physical Layer Aspects[END_REF]. The noise power is N = -96dBm.

We assume the peak bit-rate equal to 30% of the ergodic capacity of the AWGN channel; cf. Footnote 3, with the frequency bandwidth W = 5MHz and the Rayleigh fading with mean power E[|H| 2 ] = 1.

Estimations of the typical cell are performed by the simulation of 30 realizations of the Poisson model within a finite observation window, which is taken to be the disc of radius 2.63km. We first average over all BS in this window and then over the model realizations. The empirical standard deviation form the obtained averages will be presented via error-bars.

We shall study now our model using the typical and mean cell approach, assuming first the full interference model and then the weighted one.

B. Full interference

We consider first the full interference model (i.e. all BS emit the signal all the time, regardless of whether or not they serve users). Figure 2 shows the mean cell load of the typical cell E 0 [θ(0)] and the stable fraction of the network π S obtained from simulations, as well as the analytically calculated load of the mean cell θ, versus mean traffic demand per cell ρ/λ BS . We confirm that the typical cell and the mean cell models have the same load. Note that for the traffic demand up to 500kbps per cell we do not observe unstable cells in our simulation window (π S = 1).

Figure 3 shows the mean number of users per cell in the stable part of the network N 0 /π S (obtained from simulations) and the analytically calculated number of users in the mean cell N versus mean traffic demand per cell. We have two remarks. For the traffic demand smaller than 500kbps per cell (for which all the simulated cells are stable; π S = 1, cf. Figure 2), both models predict the same mean number of users per cell. Beyond this value of the traffic demand per cell the estimators of the number of users in the typical cell become not accurate due very rapidly increasing fraction of the unstable region. (Error bars on all figures represent the standard deviation in the averaging over 30 realizations of the Poisson network).

Finally, Figure 4 presents the dependence of the mean user throughput in the network on the mean traffic demand per cell obtained using the two approaches: r 0 and for the typical cell and r for the mean cell. Again, both models predict the same performance up to roughly 500kbps.

C. Weighted interference

We consider now the load-weighted interference model taking into account idling cells. We see in Figures 5,6 and 7 that the consequence of this (more realistic) assumption is that the cell loads are smaller, a larger fraction π S of the network remains stable, and the two approaches (by the typical cell and by the mean cell) predict similar values of the QoS metrics up to a larger value of the traffic demand per cell, roughly 700kbps. Note that it is in this region that the real network operates for which we present the measurements, and that its performance coincides with the performance metrics calculated using the typical and mean cell approach. More precisely, the field measurements on Figures 5,6 correspond to the same day and network zone considered on Figure 1.

Remark 14: [Measurement methodology] Measurement points on Figure 5 show the fraction of time, within a given hour, when the considered base stations were idle, averaged over the base stations, as function of the average traffic demand during this hour. Similarly, measurement points on Figure 6 show the spatial average of the mean number of users reported by the considered base stations within a given hour, as function of the average traffic demand during this hour. Finally, measurements on Figure 7 give the ratio of the total number of bits transmitted by all the base stations during Remark that Figure 7 makes evident a macroscopic relation between the traffic demand and the mean user throughput in the network zone already considered on Figure 1. This relation, we are primarily looking for in this paper, is not visible without the spatial averaging of the network measurements described in Remark 14. In order to ensure the reader that a relatively good matching between the measurements and the analytic prediction is not a coincidence we present on Figure 8 similar results for a urban zone of a different European city, where the spatial homogeneous Poissonianity of the base station locations can also be retained; cf. Figure 9. The only engineering difference of this network zone with respect to the previously considered dense urban zone is roughly four times smaller density of base stations, more precisely 1.15 stations per km 2 .

Remark 15: [Day and night hours] Let us make a final remark regarding the empirical relation between the mean user throughput on the mean traffic demand revealed on Figures 7 and8. Recall that different points on these plots correspond to different hours of some given day. In fact, the points laying below the mean curve correspond to day hours while the points laying above the mean curve correspond to night hours. This "circulation" of the measured values around the theoretical mean curve, indicated on Figure 8 and visible on both presented plots of the throughput, seems to be a more general rule, which escapes from the analysis presented in this paper and remains an open question. A possible explanation can lay in a different space-time structure of the traffic during the day and night, with the former one being much more clustered (fewer users, requesting larger volumes, generating less interference and overhead traffic).

V. CONCLUSION

In order to evaluate user's QoS metrics, in particular the mean user throughput, in large irregular multi-cellular networks, two approaches based on stochastic geometry in conjunction with queueing and information theory are developed. The typical cell approach consists in considering true spatial In fact, in large cities spatial, homogeneous "Poissonianity" of base-station locations is often satisfied "per zone" (city center, residential zone, suburbs, etc.). Moreover, log-normal shadowing further justifies Poisson assumption, cf. [START_REF] Błaszczyszyn | Using Poisson processes to model lattice cellular networks[END_REF][START_REF] Brown | Cellular performance bounds via shotgun cellular systems[END_REF]. averages of local network characteristics and thus capturing the global network performance. A simpler, approximate but fully analytic approach, called the mean cell approach, is inspired by the analysis of the typical cell. The key quantity explicitly calculated in this approach is the cell load. In analogy to the load factor of the (classical) M/G/1 processor sharing queue, it characterizes the stability condition, mean number of users and the mean user throughput. Considering some real field measurement, we validate the proposed approach by showing that it allows to predict the performance of a real network.

The present work raised open theoretical questions regarding the stability of spatially and, more difficult, space-time dependent processor sharing queues modeling the performance of individual network cells (cf Section III-D). More work is also required to figure out the problem of different performance of the network during day and night hours (cf Figure 8).

Fig. 1 .

 1 Fig. 1. Local user throughput versus local traffic demand for some zone (selected to satisfy a spatial homogeneity of the base stations) of an operational cellular network deployed in a big city in Europe. 9288 different points correspond to the measurements made by different sectors of different base stations during 24 different hours of some given day.

Fig. 2 .Fig. 3 .Fig. 4 .

 234 Fig. 2. Cell load and the stable fraction of the network versus traffic demand per cell in the full interference model.

Fig. 5 .Fig. 6 .Fig. 7 .

 567 Fig.5. Load and the stable fraction of the network versus traffic demand in the weighted interference model. Also, load estimated from real field measurements.

Fig. 8 .

 8 Fig. 8. Mean user throughput in the network versus traffic demand per surface for an urban zone of a big city in Europe. (The density of base stations is 4 times smaller than in the dense urban zone considered on Figure 7).

Fig. 9 .

 9 Fig. 9. Ripley's L-function calculated for the considered dense urban and urban network zones. (L function is the square root of the sample-based estimator of the expected number of neigbours of the typical point within a given distance, normalized by mean number of points in the disk of the same radius. Slinvyak's theorem allows to calculate the theoretical value of this function for a homogeneous Poisson process, which is L(r) = r.)In fact, in large cities spatial, homogeneous "Poissonianity" of base-station locations is often satisfied "per zone" (city center, residential zone, suburbs, etc.). Moreover, log-normal shadowing further justifies Poisson assumption, cf.[START_REF] Błaszczyszyn | Using Poisson processes to model lattice cellular networks[END_REF][START_REF] Brown | Cellular performance bounds via shotgun cellular systems[END_REF].

According to the formalism of the theory of point processes (cf e.g.[START_REF] Daley | An introduction to the theory of point processes[END_REF]), a point process is a random measure Φ = j δ X j , where δx denotes the Dirac measure at x.

Stationarity means that the distribution of the process is translation invariant, while ergodicity allows to interpret some mathematical expectations as spatial averages of some network characteristics.

The theoretical upper-bound of such function characterizing the link level performance is typically given by information theory. For example, in the case of AWGN channel, R (SINR) = W log 2 (1 + SINR) where W is the bandwidth. We shall assume that the fading is already averaged out by considering the so-called ergodic capacity. For example, in the case of flat fading,R (SINR) = W E log 2 1 + |H| 2 SINRwhere H is a random variable representing the fading and E [•] is the expectation with respect to H.

This can be achieved using various multiple access schemes, e.g. time division.

Small user movements are reflected in channel fading; cf. the remark in footnote 3.

Note that the (mean) QoS characteristics of users in this state correspond to time-averages of user characteristics.

Analysis of more sophisticated power control schemes is beyond the scope of this paper.