Heuristic Radio Resource Management for Massive MIMO in Satellite Broadband Communication Networks

The advent of technologies allowing the implementation of flexible satellite payloads is opening up great opportunities to efficiently cope with non-uniform traffic conditions. The latter is one of the “crux” of modern satellite broadband access service economic provision. In this paper we derive affordable (linear) complexity radio resource management (RRM) techniques for satellite multi-beam broadband communication networks that are closely approaching the optimum solution which has an exponential complexity with the number of users. The proposed heuristic RRM solution is applied to the most demanding Massive Multiple Input Multiple Output (M-MIMO) case. It is shown that the low-complexity payload pragmatic M-MIMO solution combined with heuristic RRM is able to efficiently cope with very different traffic conditions with a small penalty in terms of throughput compared to the much more complex ideal Minimum Mean Square Error (MMSE) M-MIMO adopting the same RRM algorithm. Instead, the pragmatic M-MIMO with heuristic RRM is outperforming MMSE M-MIMO exploiting simple RRM.


I. INTRODUCTION
The need for higher throughput, hence large frequency reuse, in point-to-point satellite communications has been triggering the exploitation of payloads with a large number of beams and the development of associated technologies. While the increase of the number of satellite beams is beneficial in terms of achievable frequency reuse and antenna gain, its practical adoption creates a number of system-level issues that have been triggering researchers' interest. The main challenge is related to the non uniform-traffic nature seen by a satellite due to its large coverage area. This is because the satellite broadband access services are playing a complementary role to terrestrial networks. Hence, satellite's service penetration is typically related to under-served regions such as low-densely populated areas, maritime and aviation routes over non populated areas (e.g. oceanic routes). To cope with this spatially uneven and time variant traffic distribution over the satellite coverage area, a number of technologies to support flexible payload resource allocation have been devised and are being continuously improved. Among the most important ones, we mention digital on-board processors, beam hopping, adaptive coding and modulation, flexible power allocation per beam and active antennas [1]. Satellites equipped with flexible payloads and resource allocation, thanks to their high altitude and large coverage area, have the potential to mitigate the demand-supply space-time mismatch. On the contrary, terrestrial networks deployment is sized for a given traffic pattern (cell size and density) with limited short/medium term flexibility to cover demand evolution.
As the flexible multi-beam payload technologies were becoming available at affordable cost, the problem on how to manage the intrinsic space segment flexibility aroused. As a consequence, Radio Resource Management (RRM) for multi-beam satellite networks (with or without precoding interference mitigation techniques) became an active field of research in the last fifteen years as the large relevant literature testifies [2]- [22].
Previous work on RRM for multi-beam satellite networks has been assuming, with few exceptions, a conventional geostationary orbit (GSO) with single feed per beam payload architecture. Therefore the above references are not relevant to the case of M-MIMO satellite networks.
An overview of RRM research for terrestrial networks Multi-User MIMO is reported in [16]. In this paper it is mentioned that RRM for M-MIMO is an open research field. Regarding satellite systems, in [23], a M-MIMO for LEO satellites simplified approach is investigated. In particular, to make the precoding easier to implement at physical layer level, it is proposed to adopt a statistical instead of instantaneous channel information. In addition, the authors propose a space angle based grouping to schedule together users into different groups. The proposed approach is shown to be asymptotically optimum for a sufficiently large number of satellite antenna feeds and user groups.
In this paper we extend the work initiated in [24], looking at an heuristic computational efficient solution for satellite multi-beam pragmatic M-MIMO approximating the Mixed Integer Quadratic Programming (MIQP) approach described in that work. It is shown that the proposed heuristic RRM (H-RRM) algorithm has performance very close to the MIQP with a complexity linearly (instead of exponentially) dependent on the number of users. The H-RRM solution optimally organising the active users' population in the system's resource colours (slices in time, frequency or code) can be further enhanced by a further beam selection optimization step based on the so called "virtual distancing" technique. The heuristic RRM solution proposed can be easily implemented using the current DVB-S2X standard and perfectly matches the Multi-Beam (MB) M-MIMO simplified payload architecture proposed in [24]. Although the proposed RRM technique is particularly suited for the pragmatic M-MIMO case it can also be used in more conventional multi-beam satellite systems exploiting frequency colouring schemes.
The two new "virtual distancing" approaches allow to determine beam steering positions and relevant precoding/beamforming vectors based on the knowledge of the contemporary active users with some improvement of the overall system throughput. Although they are described in relation to the pragmatic-MIMO and H-RRM approaches, they have broader applicability to multibeam flexible payloads.
Both the H-RRM and the virtual distancing techniques only require approximated (e.g. few km accuracy) information of the user geographical location, the antenna beam pattern and imply a low implementation complexity.
Through extensive computer simulations the H-RRM techniques have been applied to the pragmatic M-MIMO multibeam satellite system in the presence of many non-uniform and uniform (best case) traffic scenarios. For all traffic distribution cases investigated, the pragmatic M-MIMO throughput performance are shown to be only slightly below to ideal Minimum Mean Square Error (MMSE) M-MIMO precoding with ideal channel estimation and adopting the same H-RRM scheme. Instead, when using a simple RRM scheme for MMSE M-MIMO, the performance of the latter are well below the pragmatic M-MIMO using H-RRM. Simulation results are showing a truly remarkable H-RRM algorithm robustness to very different traffic scenarios and allowed to understand the system throughput performance dependency on key system parameters. It is remarked that, as discussed in [24], MMSE M-MIMO represents a purely theoretical upper bound as it can not be implemented in practice without major impairments. Instead, the proposed H-RRM combined with the pragmatic M-MIMO approach represents a readily feasible solution.
The paper is organized as follows:in Sec. II the M-MIMO system is described in terms of forward link, physical layer and traffic models. Sec. III outlines the proposed H-RRM algorithms approximating the MIQP solution and the virtual distancing techniques. Sec. IV provides extensive simulation results about the H-RRM and virtual distancing algorithms under different traffic scenarios. Finally, Sec. V provides the summary of the paper and associated conclusions.

II. SYSTEM MODEL A. FORWARD LINK MODEL
In the following we will assume that a population of N users has to be served by the satellite forward-link. Assuming that a precoding/beamforming technique [24] is applied to the full population of users 1 can be introduced such that its entry element s(i, j) represents the power transmitted for the j-th user and received by the i-th user, normalised to the thermal noise power. Without loss of generality it is assumed that all the users have identical terminal characteristics and experience same receive noise power. We will use the following notation: vectors will be represented in bold lower-case and matrices in bold capital letters. Superscripts T and H will indicate transpose and complex conjugate transpose matrix operators, respectively.
Signal-to-noise ratio (SNR), interference-to-noise ratio (INR) and signal-to-noise plus interference ratio (SNIR) for the i-th user can be derived from relevant entries of the power transfer matrix matrix S, respectively .
The diagonal entries of the matrix S thus represent the useful part of the transmission corresponding to the signal-to-noiseratio. All the off diagonal elements are unwanted contributions that should be minimised. Taking the elements of the i-th row, with the exclusion of the diagonal one, we have the aggregate interference-to-noise ratio experienced by the i-th user. While taking the i-th column, the off-diagonal elements represent the power of the patterns for the i-th user leaking in undesired directions normalized to the noise power. The matrix S is generally non-symmetric.

B. PHYSICAL LAYER MODEL
Following [24], we assume exploiting full frequency reuse among the currently active users' beams. These are the ones simultaneously activated by the RRM algorithm and corresponding to a subset of the family of beams that the payload active antenna can generate. As only a subset of the potential beams generated by the active antenna are active at the same time, a beam colouring scheme has to be introduced in the physical layer. The proposed colouring scheme is based on organising the different sub-set of beams transmission in Time Division Multiplexing (TDM) with a super-frame (SF) having time duration T SF . The preference for TDM instead of Frequency Division Multiplexing (FDM) is related to its easier payload implementation. In fact,TDM although fully dual to FDM in terms of performance, allows to avoid the flexible payload channelization required by the FDM approach. TDM can also be combined with an FDM approach to reduce the TDM carrier bandwidth. In this case the RRM should be preferably applied independently to each TDM sub-band as the payload may have a frequency dependent beam forming characteristic over a very wide bandwidth. An alternative colouring scheme is represented by Code Division Multiplexing (CDM) which is also allowing full frequency reuse. However, CDM is less used in satellite networks and standards, therefore will not be considered in the following. As shown in Fig. 1, the SF is divided in N FR frames (FR) having duration T FR . Each FR represents a colour to be used by the RRM to place users at the same time over the satellite coverage region. Each FR is divided in N SL slots (SL) having duration T SL . To each slot the RRM assigns a specific Adaptive Coding and Modulation (ACM) physical layer configuration called MODCOD in line with the DVB-S2X standard nomenclature [27]. The FR partitioning in SL ensures a high level of FR filling. This is because, as we will see, each frame is associated to a specific satellite beam. Allowing different MODCODs inside each FR allows to serve different users within the current beam experiencing different signal-to-noise plus interference ratios. Each SL is divided in N SY symbols each having duration T s , being R s the TDM carrier transmission baud rate, with R s = 1/T s . The following relations hold For the SL s in frame number f the number of bits transmit- being η ACM (s, f ) the spectral efficiency corresponding to the MODCOD selected for the slot s in frame f . The total number of bits transmitted in the FR number f is simply given by Similarly for the SF the total number of bits transmitted is The corresponding SF average bit rate is The number of time slots in a frame depends on the SNIR range to be covered by the ACM MODCODs and the number of MODCODs supported. Typically each MODCOD is carrying packets belonging to different users sharing similar received SNIR values and served by the same satellite beam. Adopting a large number of MODCODs reduces the ACM SNIR granularity. However,this approach is likely to diminish the filling rate of each MODCOD. Consequently, a trade-off is required to optimize the overall physical layer efficiency. In the following, we will look at the RRM in terms of colour assignment per FR with ACM assuming there is enough traffic request to achieve a full SL packets filling rate 2 . As it will be clear in the following H-RRM technique detailed description, the proposed approach does not require to adapt the super-frame format to the actual traffic distribution. This simplifies the beam-hopping functioning and associated signalling implementation.
An important aspect is related to the ACM implementation. The pragmatic M-MIMO approach described in [24] implies the exploitation of beam hopping (BH) to activate only a subset of the satellite beams at any time. As it will be shown in Sec. III, the proposed RRM algorithm will use the BH payload capability in a dynamic way to serve the traffic request. This means that the subset of active beams will be changing on a frame-by-frame basis but repeating on superframe-by-superframe basis at least for a reasonable amount of time 3 . To adapt the ACM physical layer configuration [25], the SNIR evaluation shall be restricted to the useful user frame. This is to avoid unwanted fluctuations of the ACM SNIR estimates generated by other frames dedicated to users located in nearby active beams of the SF.

C. TRAFFIC MODEL
To test the H-RRM algorithm performance under various traffic conditions we have been developing a flexible spatially non uniform random traffic model. The model is assuming that the traffic is present in N TS circular spots. The generic traffic spot l is characterized by a normalized radius r TS (l) and center location normalized offset distance with respect to the satellite nadir c TS (l) with orientation angle φ TS (l). Each traffic spot has a probability of occurrence p TS (l), with NTS l=1 p TS (l) = 1. Inside each spot where the traffic can generated with Gaussian bi-dimensional truncated distribution within the circular spot area characterized by the standard deviations σ XX , σ YY and σ XY or with an uniform distribution. The simulated traffic configurations are summarized in Tables 1-2 for the uneven and even spot traffic probabilities respectively.
Case 1 corresponds to a single traffic area S T covering the full satellite coverage area S A with uniform traffic distribution (or Poisson Point Process see Appendix A). In case 2 we have 5 traffic spots which in sub-cases a and b are overlapping among them so that S T = S A . Instead, for sub-cases c, d, e, f we have S T < S A . Case 3 corresponds to a 2 traffic spots with different geometry and traffic distributions in the sub-cases a-f. In all cases offset beams are characterised by c TS (l) = 1/3.
The simulated PDFs for the Table 1 geometry configurations are shown in Fig. 2, 3, 4, and Table 2

III. RADIO RESOURCE MANAGEMENT ALGORITHMS
Considering that the user distribution is an uncontrollable input with non-uniform statistical distribution and possible clustering due to hot spots in areas of high user density we described in [24] a method to assign the users to different radio resources (time slots and/or frequency carriers) in a manner that maximise the overall throughput. We recall that the generic element s(i, j) of the power transfer matrix S introduced in [24] describes the power transmitted for the user j and received by the user i. Here for convenience the matrix S has dimension (N ×N ) where N represents the total number of active users in all colours or time slices. Instead N U (c) represents the number of users allocated by the RRM to the colour (or slice) number c. To be remarked that, as it will be shown in Sec. IV-C4, N U (c) is a integer random variable when considering different SFs. Disregarding for simplicity of notation the N U dependency on the SF index, for each SF the following relation holds Analysing the power transfer matrix S, it can be easily recognised that the off-diagonal coefficients s(i, j) play a detrimental effect in increasing the interference experienced by the i-th user. Assigning the users to different time/frequency slots nulls their reciprocal interference. This effect can be modelled with a generalised co-channel matrix X for which a unit entry identifies the assignment of the same radio resource to the users with relevant row and column indexes as x(i, j) = 1 if user i and user j share the same colour, 0 otherwise.
The matrix X is symmetric with unit diagonal elements and has the same dimensions of the power transfer matrix S. They can be combined by mean of the Hadamard product in a resulting co-channel power transfer matrix S CC which takes into account the radio resource assignments. Following [24] the matrix S CC can be computed as The SNR, INR and SNIR equations for the i-th user (i = 1, · · · , N ) are modified accordingly as Given a power transfer matrix S the resource management problem is cast in a colour assignment problem which is formulated as a 0-1 integer programming instance on the colouring matrix binary entries. Defining C as the preassigned number of colours, the binary colouring matrix C of size (N × C) can be introduced Considering that the i-th user can be assigned only to one colour, the rows of C must satisfy the following linear constraint corresponding to one colour per user Additionally, summing all entries The co-channel matrix X can be obtained as the matrix product of C with its transpose C T Considering that the central target of the colour assignment is the reduction of the mutual interference we focus our attention on the off-diagonal elements of the power transfer matrix introducing its diagonally null part Q which we briefly term the interference matrix where the diag {·} operator is used accordingly to MATLAB notation (i.e. returning a vector with the diagonal elements of a square matrix, and returning a diagonal square matrix if applied to a vector). With the nulling of the diagonal elements the interference-to-noise ratio for the i-th user (INR(i)) summation can be extended to all the indexes 4 x(i, j)q(i, j).
The sum co-channel interference, evaluated as the total interference-to-noise ratio INR T , can be assumed as a figure of merit of the overall performance and satisfies the following definition In matrix form we can write INR T in different equivalent forms (refer to [24])  (25) While several techniques have been proposed for solving in exact or sub-optimal manner MIQP problems, an efficient heuristic approach to a strictly Integer Quadratic Programming (IQP) problem has been proposed in [26] for a beam colouring problem. A solution to that problem requires that adjacent beams do not share the same colour, which can be recognised as a vertex coloring problem of graph theory. The search of a satisfactory solution is addressed in [26] with a heuristic iterative colouring of the beams which satisfy the adjacency constraints. By reinterpreting the binary matrix representing the adjacency conditions between beams as a binary interference matrix, and casting the colour compatibility condition in a quadratic form, we can translate the beam colouring problem in a sum-interference minimization problem for which we seek a null sum-interference solution. Generalizing the binary adjacency matrix of the beam-colouring IQP problem [26] with the real positive interference matrix with null diagonal of the sum-interference MIQP [24], and introducing an auxiliary interference matrix which provides a quantitative selection criteria at each iteration we can extend the heuristic search algorithm to the MIQP problem at hand. The proposed generalization is discussed in the following section.

A. HEURISTIC RRM
We introduce the auxiliary co-channel interference matrix A of size (N × C) such that a(i, k) represents the aggregated co-channel interference on user i due to the co-frequency interference of user in colour k. The element a(i, k) can be written as the product of the i-th row vector q R i of the matrix Q times the k-th It is worth noting that the matrix A contains an excess of information as each of the entries of the i-th row reports the interference on the i-th user as if it would be associated to the k-th colour, with a number of entries equal to the number of colours k = 1 . . . C.
Left multiplying the matrix A by C T we obtain a square matrix of size (C × C) with entry elements  which sum-up the interference of user assigned to colour l as they would have experienced interference from users in colour k. Obviously this quantities make sense only when the indexes l and k coincide. It this respect the diagonal elements C T A (k, k) represent the sum interference for colour k.
The sum-interference INR T can be finally obtained summing all the diagonal elements The richness of information offered by the matrix A can be exploited in a constructive way to assign colour to users in an heuristic sequential way. Assuming that at the n-th iteration a number of users has been assigned to a colour such that the temporary colouring matrix for the non-assigned users, the information provided by A (n) = Q C (n) allows to select the one (row) with worst interference and for it the least interfering colour (within the row). Using the concepts above, the pseudo-code for the first colouring iteration is described in Algorithm 1. In doing so we give priority to the most interfered users in picking up the least interfering colour. This is the key strategy beyond the proposed heuristic approach that the matrix formulation allows to be easily implement in software. As per pseudocode, the first user selected and assigned to colour 1 is the one showing the highest interference from any other user.
In the evaluation of A (n) = Q C (n) it can be observed that we need to evaluate it only for the non-assigned users, user_flag(i) = 0, considering that the interference comes only from the assigned users user_flag(i) = 1. This allow to dynamically optimize the size of the matrices to the smallest dimension with relevant multiplicative complexity saving.
For the described algorithm, the combinatorial complexity remains strictly linear O(N ).

B. FAST HEURISTIC RRM
The heuristic colouring algorithm has been described and implemented in Algorithm 1 with a matrix notation pivoted in the iterative evaluation of the user interference by mean of the matrix multiplication A (n) = Q C (n) . This matrix multiplication dominates the computational complexity of the algorithm 5 with a multiplicative complexity of order A number of improvements on the data structure and data manipulations can be implemented to increase the efficiency of the evaluation of the matrix product: • Considering the binary nature of the matrix C (n) , the matrix multiplication may be interpreted as a selective sum with the indexes of the elements to sum defined by the entries of the column vectors of the matrix C (n) . This leads to computational improvements due to the avoidance of multiplications. • Exploiting the sparse nature of the interference matrix Q, which is due to the decaying nature of the radiation patterns, a data structure based on lists of non-null elements, instead of a matrix structure, could be more effective. This approach is already well established in graph theory literature which addresses similar combinatorial computational problems (refer for example to Chapter 22 of [28]). The matrix Q can be interpreted as a weighted adjacency matrix which can be efficiently represented by its adjacency list. • Alternatively, a similar improvement in the matrix multiplication can be obtained with sparse matrix representation techniques [29], [30]. The matrix Q can be reordered in a band matrix Q (i.e. a matrix whose only nonzero elements lie on diagonal bands above and/or below the main diagonal) of limited bandwidth (i.e.  the number of nonzero diagonals below and above the main diagonal) by means of matrix bandwidth reduction techniques, such as the Cuthill-McKee algorithm [31], or its reverse version [29]. The reordering of Q in Q has a limited combinatorial complexity O(N 2 ) [32] and needs to be performed only at the beginning of the heuristic user colouring routine. The effect of the application of the reverse Cuthilla-McKee algorithm the for a population of 200 users generated according to a random uniform distribution in the coverage area is shown in Figure 7. The graph on the left side of the figure shows the distribution of the non-negligible elements (i.e. below a given threshold) of the approximated power transfer matrix; on the right side the reordered bandwidth limited matrix. The matrix multiplication A (n) = Q C (n) can be performed limiting the running indexes of the entrywise row-column multiplication to non-null elements within the identified bandwidth of Q . The elements above address the possible data structures that can be exploited to improve the efficiency of the matrix multiplications A = Q C. Nevertheless, a dramatic reduction of the computational complexity can be obtained observing in further details the operations involved in the transition between two iterations. Lets focus our attention on the result of iteration (n) which is used as starting point for the matrix multiplication at the following iteration (n + 1) The binary colouring matrix C (n+1) is the result of an update of the colouring matrix C (n) , which is the input of iteration (n) with the assignment of colour k (n) to user i (n) during iteration (n). It corresponds to setting to one the element of C (n+1) with row and column indexes i (n) , k (n) The evaluation of the new auxiliary co-channel interference matrix A (n+1) can be thus reduced to an update of the matrix Considering the structure of ∆C (n) , the matrix multiplication ∆A (n) = Q∆C (n) corresponds to selecting the i (n) -th column of Q, q C i (n) , and placing it as k (n) column of the delta This result can be interpreted as the fact that with the assignment of the colour k (n) to the user i (n) , its interference (i.e. q C i (n) ) must be added to the relevant colour column (i.e. k (n) ) of the auxiliary co-channel interference matrix. This allows to reduce the evaluation of A (n+1) to an update of the matrix at previous iteration, A (n) , adding the interference of the newly assigned user to the column of the selected colour, ∆A (n) , with a computational complexity of N real additions in the worst case. In this respect, a representation of the matrix Q as a re-ordered band-limited matrix or with an adjacency list would allow to further reduce the number of additions either to the number of non-null entries per column or to the vertical band of the matrix, respectively.
A fast version of Algorithm 1 which exploits the iterative update of the auxiliary co-channel interference matrix is reported in Algorithm 2. It is worth noting that for economy of operations, also the completion condition of the while loop is expressed with an iteratively updated counter in both algorithms.

C. APPROXIMATED POWER TRANSFER MATRIX
Depending on the beamforming/precoding approach the matrix Q, which collects the mutual interference terms, can be evaluated in an approximated or detailed manner. An a priori evaluation of the interference matrix Q based on the mutual distance is possible and is proposed as a an effective solution to the resource management problem.
Considering that main contributions for the interference are due to main beam, a Gaussian shaped beam g 0 (u, v), VOLUME 4, 2016 Algorithm 1 Heuristic User Colouring Inputs: S , Power transfer matrix s(i, j), i, j = 1 . . . N C , Number of colours N , Total number of users (i.e. in all colours) Can be derived from the size of S

Outputs:
C , Colouring matrix Find user with worst aggregated interference Removes user i (n) from the list end while of beam-width comparable with the expected antenna beamwidth, can be used as reference prototype beam The parameter σ can be adjusted to to match the main beam of the antenna under assessment, which in [24] is a square array of size D A . For simplicity we use the expression of the equivalent continuous square aperture g S (u, v) Equating the two expressions at the point u = 1 Can be derived from the size of S

Outputs:
C , Colouring matrix Counter of assigned users Find user with worst aggregated interference Removes user i (n) from the list A(:, k (n) ) = A(:, k (n) ) + Q(:, i (n) ) user_counter ← user_counter + 1 end while we obtain and finally The Gaussian beam pattern approximation has a 4 dB beamwidth corresponding to the array of size D A /λ. A comparison of the two pattern under the main-beam matching condition is reported in Fig. 8.
The normalised interference term q(i, j) can be thus evaluated as In general, it is recommended to use the full antenna beam pattern as it provides the most accurate RRM solution with an affordable complexity. Nevertheless, in the rest of the paper the Gaussian beam pattern approximation will be used only for the RRM allocation of users to colours, while the capacity evaluations will be based on the full analytical antenna model detailed in [24].

D. HEURISTIC RRM ALGORITHM EXEMPLARY RESULTS
In this section we report some exemplary application of the heuristic user colouring algorithm. It is worth noting that Algorithm 1 and 2, respectively, have the same combinatorial decision structure and apart negligible numerical error propagation effects they shall provide the same outputs for the same inputs. The heuristic user colouring Algorithm 1 and 2 have been applied to a uniformly random distribution of M = 600 users within the (u, v) circle defined by the Earth as seen from a Geostationary satellite. The random user distribution is shown in Fig. 9-a. Application of heuristic user colouring (Algorithm 1 and 2) leads to the colouring shown in Fig.  9-b, with the individual colours slices shown in Fig. 10. The circles indicates 4 dB beam-width contour plots for an uniformly illuminated square antenna of size D A = 1.2 m at f = 20 GHz. It can be observed that the populations of users per colours have an increased average separation.

E. HEURISTIC RRM ALGORITHM STOCHASTIC ANALYSIS
Insightful information on the performance of the H-RRM Algorithm can be gained analysing it from a stochastic geometry perspective. In mathematics, stochastic geometry is the discipline studying spatial patterns [33], [34] and it has been extensively applied to modelling of terrestrial wireless networks performance [35], [36].
A spatial point pattern can be defined as a set of point locations within a Region of Interest (ROI). Each point, also indicated as event, is identified by its coordinates. For our purposes we can consider the user distribution in the u, v plane as a point pattern. Application of the H-RRM Algorithm results in a number of slices, one per colour, and we can analyse the geometrical properties of the slices to obtain performance trends.
The benchmark model for spatial point patterns is the Poisson Point Process (PPP) (see Appendix A) which exhibits Complete Spatial Randomness (CSR) [33]. In this model the points follow a Poisson Process over the study region which results in a uniform probability density function of intensity ρ, being ρ the mean number of points per unit area.
As observed in [24], optimal throughput performance can be obtained when the user distribution satisfies a nearest neighbour minimum distance (referred to as Poisson Disk Process (PDP) distribution, see Appendix B) roughly equal to the antenna pattern resolution. It is so arguable that we may focus our attention to first-order statistics of the slice point pattern, in particular on the statistical distribution of nearest-neighbor event-event distance.
The nearest-neighbour distance distribution functionĜ(r), known as Diggle's G-function [37], describes the cumulative distribution of the distance of a random point to the nearest point (event-to-event) in the point pattern. A comparison of the empirical distribution function of the point pattern under analysis with the theoretical distribution of the PPP with the same density of points (for details on the normalization rescaling refer to Appendix A) offers a summary of effects that the RRM induces on the distancing of the points of a same slice. Figure 11 shows the dependency of the empirical distributions of nearest-neighbour distances for different sizes of populations per slice (N/C = 250, 500 and 1000) and different number of slices (colours). The re-scaling of the empirical distributions to a normalized density evidences that the effects of the H-RRM are substantially independent on the size of the population and depends mainly on the number of slices (colours).
The asymptotic effects of the number of slices (colours) on the distributions of nearest -neighbour distances can be appreciated in Fig. 12 where a saturation effect can be observed. The figure reports also a best fitting of the empirical distributions with the cumulative distributions of Weibull three-parameters PDFs. The best fitting parameters have been obtained by Maximum-Likelihood Estimation [40].
Additionally, the following aspects are worth of note: • The Algorithms 1 and 2 have been described as a applied to the full user population. Tailored versions of them can be used to assign colour to new users entering in the active population, or to reassign resources freed by users becoming inactive.  • It can be noted from the exemplary results that the minimum distance (i.e. disks not intersecting) cannot be guaranteed in all the cases when a given user population and number of colours are assigned. This can be justified by the observation that applying a moving spatial filter of the size of the antenna disk, spatial locations counting a number of users within said disk higher than the number of colours will not be resolvable. • Equalization of the number of users assigned per colours is not a direct target of the algorithm. Nevertheless, recalling that INR T = tr C T QC , and that the diagonal entries of the matrix C T QC of size (C × C) represent the sum-interference per colour, it is observed that the heuristic algorithm tends to equalize the value of these diagonal entries and, as secondary effect, obtains a balanced number of users with a small deviation. This conjecture found confirmation by the simulation findings reported in Sec. IV-C4.

F. VIRTUAL DISTANCING TECHNIQUES
In order to enhance the previously described H-RRM algorithms, some area of improvement has been identified in tackling the users' beam assignment worst-cases. The virtual distancing techniques described in the following two sub-paragraphs are based on the observation that once the colours are assigned, the beam positions determine the useful gains and unwanted interference for all users in the resource slot (colour) under consideration. An optimization of the beam positions which guarantees a useful gain for the wanted user while reducing the interference to the users in the same slot may obtain an improvement in the overall throughput.
Considering that the main-lobe interference is mainly dependent on the distance of the beam position to the inter-  neighbour users. In this manner the distance between the users is "virtually" increased.
Two different approaches will be described in the following sub-paragraphs and their performance analysed in Sec. IV: • Virtual Distancing (VDI) -based on the optimization of the beam positions as continuous bi-dimensional variables; • Quantized Virtual Distancing (QVD) -where the possible beam positions are discrete points of a bidimensional lattice which has payload simplification effects [24].
Both approaches benefit from a form of uplink-downlink duality [41] that allows to optimize individually the beam positions, so reducing the overall complexity of the problem.
This can be made more explicit observing that following (24) the sum-interference INR T can be expanded in the sum of the elements of the main diagonal of the matrix C T QC The square matrices Q k = c C k T Q c C k are sub-matrices 6 of Q for colour k. They are obtained retaining the rows and columns of Q corresponding to the active users for colour k. Left and right multiplication respectively by 1 T and 1 (where vector 1 has all unit entries), returns the sum of all the matrix entries of Q k (being the matrix entries positive, it corresponds also to the L 1 matrix norm).
The associative property of matrix multiplication, 1 T Q k 1 = 1 T (Q k 1) = (1 T Q k )1 , allows to interpret each colour contribution to INR T either a sum on downlink interference Q k 1 or uplink interference 1 T Q k . To better understand the difference we consider for simplicity the case of a single colour where Q k can be read as Q. The downlink and uplink interpretations follow: The column vector Q1 has row entries, which correspond to the sum of the interference of users j on user i sharing the same colour. So each row element of the column vector Q1 can be interpreted as a downlink interference per user; all the users of colour k contribute to the interference experienced by user i.

• UPLINK INTERPRETATION
The row vector 1 T Q has column entries, Each column element of the row vector 1 T Q can be interpreted as the aggregated interference generated by user j on all the other users. Similarly to the uplink case, the user j is interfered by all the users sharing the same colour.
Representing that the i-th user position in the u, v-plane with u i and assuming that the j-th beam is pointing at s j , where where s x j and s y j are the u and v coordinates of the beam pointing direction, respectively, the interference generated by user j on user i can be considered mainly determined by the 6 A matrix obtained by deleting some of the rows and/or columns of a matrix is said to be a sub-matrix of a given matrix. VOLUME 4, 2016 spatial response of the j-th beam array factor where g A is the nadir-pointing array factor.
Exploiting the uplink interpretation, the minimization of the sum-interference can be translated in the independent optimization of the beam positions s j such that, for each beam, the aggregated interference is minimized. It should be recalled that beam s j serves user u j . To limit the gain loss due to the beam displacement from the optimal gain position, the optimization search space must be constrained to a maximum gain loss ∆G max Considering that the un-steered array factor g A (·) has in 0) a maximum and stationary point, the constraint on the maximum gain loss can be translated in a constraint on the a maximum displacement radius of the beam pointing s j Algorithms for continuous or quantized optimization of the beam positions are described in the following paragraphs.
To simplify the lexicographic use of the indexes involved in the description of the optimization algorithms, the running index of the beam under optimization will be denoted as i and the reader is warned of the possible confusion that the index exchange may induce.

1) Virtual Distancing (VDI) Algorithm
In [24] the authors reported results of sensitivity analyses of the system capacity on the normalised minimum distance ρ min between users 7 . The parameter ρ min is defined as the Euclidean distance between two points in the u, v plane, normalised to the approximated array beam-width λ/D A (as u, v diameter). These results suggest that the optimization of the i-th beam position s i can be translated in the analogue one of optimal distancing beam s i from the position of the co-channel users u j while limiting the gain loss for the served user u i . The Virtual Distancing (VDI) Algorithm described hereafter exploits analogies with molecular dynamics. At initialization, the beam positions correspond to the user positions (refer to Fig. 13-a. s Beam and user positions are considered as a system of particles, while beams are free to move, users remain fixed. During the optimization the beam center s i moves according to repulsive forces between the beam center and the cochannel user positions u i , i = j (refer to Fig. 13-b ). The 7 More recent simulation results showed an optimum normalized distance value of ρ min = 1.0.
force on s i due to the interaction with u j is indicated as F i,j . It has a direction corresponding to the unitary vectorr i,j , and amplitude F i,j , where The force magnitude is modelled accordingly to a modified Hooke's elastic law where: • d min is the target minimum separation between the beam and the users. It can be expressed as the product of the target normalised minimum distance between users times the approximated array beam-width • α is the spring constant. For our needs it can be assumed to be unitary α = 1. • βd min is the spring rest length. The constant β ≥ 1 is introduced to have a non negligible force for r i,j approaching d min . The results described below have been obtained with β = 1.1. The elastic force (52) is truncated to the repulsive side (i.e. F i,j ≥ 0 for r i,j ≤ d min ) so that a beam at a distance greater than d min from a user is not subject to any force. The total force acting on beam s i can be obtained by linear superposition of the individual interaction forces The position at iteration n, s (n) i , is updated by integration of the motion equation in the simplifying hypothesis of zero velocity (high viscosity condition) being ∆s i . From the mechanical point of view the constant γ covers the effects due to the particle mass and time interval. For our aims, it can be used as a dumping constant which can be tuned and/or made adaptive to control the convergence of the optimization.
Once the position difference ∆s The algorithm terminates when the displacement is below a defined tolerance. Figure 13-c shows how starting from the user positions of Fig. 13-a, the VDI algorithm converges to beam positions which fulfill the distancing condition. The final beam layout, with representative beam-widths is depicted in Fig. 13-d.
The pseudo-code for the VDI is described in Algorithm 3. Exemplary simulation results are reported in Figures 14. Figure 14-a shows maximum gain beams for a uniform distribution of users (N U = 100) on the visible earth as seen from a geostationary satellite with a square antenna of width D A = 1.2 m working at 20 GHz. The outcome of the VDI algorithm with the same user distribution of Fig. 14-a is reported in Fig. 14 Theoretical and implementation aspects of the generation of a periodic lattice of beams from a planar array are addressed in [43] and [44], respectively. The beam to user assignment can be performed according to different criteria. If maximum gain is the selected criteria [24], the bi-dimensional quantizer is defined by the Voronoi partition of the space into non overlapping cells (refer to Fig.

Algorithm 3 Virtual Distancing (VDI)
Inputs: u j , User coordinates in u, v plane j = 1 . . . N U N U , Number of users (i.e. in the colour under analysis) i , index of the the moving beam d min , target minimum distance between beam center and interfered users rmax , maximum distance between beam center to served user d tol , Termination tolerance N iter , Maximum number of iterations β, Constant to render the elastic force non-negligible close to rest (← 1.1)

Outputs:
s i , Beam center coordinates in u, v plane Initialize: s Determines user to beam versors Total force evaluation γ = 10(1 + n/N iter ) −1 Empirical adaptive damping coefficient 15-a). The i-th user, with position u i is assigned to the closest beam position s i In view of the sum-interference minimization criteria and based on the uplink interpretation, the selection of optimal beam position s i on the discrete beam lattice Λ(S) which minimize the interference to co-channel users u j , j = i, while guaranteeing a limited gain loss impact for the served user u i . The quantity to be minimized is then the interference to useful signal power ratio for user i that can be described (58) The maximum distance constraint, |s i − u i | ≤ r max drastically limit the search space in Λ(S) to the lattice points within the circle of radius r max and center u i (blue points in Fig.  15-b). For such limited number of candidates, the evaluation of the array factor for all the co-channel user locations is an affordable computational complexity.
Similarly to the procedure described in [24], the lattice of beam positions can be generated according to an hexagonal lattice base matrix S MB where the beam-to-beam interdistance is controlled by a normalised beam spacing S n MB referred to the array normalised beam-width λ/D A The pre-defined set of beams centers {s k } can be evaluated as the intersection of lattice of points Λ(S MB ) and the ROI in the u, v plane. For example, if the ROI is the disk of radius sin[ϑ max (α min )]), the set of beams centers is For the i-th user the candidates beams can be found in the subset of {s k } which fulfill the maximum distance constraint {|s k − u i | ≤ r max }. For each of these candidate beam centers the sum-interference can be evaluated and the minimum found. The pseudocode of the resulting Quantized Virtual Distancing (QVD) is described in Algorithm 4.

A. SIMULATOR DESCRIPTION
To analyze the performance of the proposed H-RRM algorithm in a multibeam satellite system context, we have been evolving the M-MIMO satellite multibeam simulator described in detail in [24] to include the required additional functionalities. More specifically. the following features have been introduced: • TDM frames 8 in the physical layer.
• H-RRM algorithms in alternative to the PDP minimum distance optimum approach. • Sec. II-C non-uniform traffic generation in addition to Appendix A PPP as input to the heuristic MIQP RRM algorithm. As explained in Sec. III, we do not have anymore a fixed amount of users allocated to each slice N U as this parameter is fluctuating depending on the specific non uniform traffic realization and consequent RRM users allocation per frame. In the following we will indicate with N U the average number of users/frame.
Although the simulator is supporting all the H-RRM algorithms described before, in the following we will not report VDI results as the performance were found to be similar or worse than the QVD ones which is also much simpler to implement in a practical satellite payload.

B. SYSTEM LEVEL ASSUMPTIONS
In line with [24] for deriving numerical results we adopted a Geostationary Satellite (GSO) configuration with the key system parameters reported in Table 3.

1) Heuristic versus MIQP RRM
As it was mentioned before, the MIQP RRM implementation has a very high complexity, hence it has a limited practical use. However, to show the goodness of the proposed H-RRM approach we have been running a small scale simulation with N = 100 users on the system model previously described and with an array size D A = 1.2 m. For the MIQP implementation we resorted to the Gurobi optimizer. The users' spectral efficiency CDF simulation results assuming Shannon physical layer performance are shown in Fig. 16. Figure 16-a shows the wide performance difference between not using RRM and the almost identical spectral efficiency obtained exploiting MIQP and H-RRM. The negligible MIQP performance advantage can be observed in the zoomed version of the same plot shown in Fig. 16-b. In the following only H-RRM solutions will be considered as the MIQP computational complexity does not allow to simulate large number of users.
Determines sum interference for user i assigned to beam n Determines useful power for user j assigned to beam n sum_I_over_C(n, i) ← sum_I(n, i)/C(n, i) Determines sum interference over useful power for user i assigned to beam n end for Find the index of the minimum of sum_I_over_C(n, i) beam_idx(min_idx) is the optimum beam index i to be assigned to user i end for

2) Impact of the Monte Carlo randomization
As it was explained in the original pragmatic M-MIMO work [24], to achieve meaningful performance results a number of Monte Carlo traffic generation trials N trials are required. Experimentally we found that 5 trials are good enough to achieve a good throughput accuracy. As an example, for the case 2f of Table 1, the MB throughput for N trials = 5 amounts to 138.85 Gbps. By increasing N trials to 10 the simulated throughput amounts to 138.91 corresponding to a negligible impact on the result. When reducing the number of trials to 1 the throughput is also similar i.e. 138.17 Gbps. These results have been confirmed for other traffic configurations showing a very minimum impact on the number of trials on the overall throughput estimation. For this reason, when not stated otherwise, for the following simulations we adopted the value N trials = 5.

3) Impact of the number of TDM frames
As we explained in Sect. II-B, the number of frames in the superframe allow to reorganise the users served at a given time to approximate the MIQP distance distribution. VOLUME  The MIQP algorithm is representing a lower bound to the optimum PDP minimum distance users' distribution used in the previous work [24]. Intuitively the greater is the number of frames N FR , the easier the H-RRM algorithm task will be. However, by increasing N FR we increase the superframe duration, hence the system latency. In Fig. 17 we have been reporting the simulated throughput as a function of the number of frames for the traffic cases 2f and 3f of Table 1 and 1.2 m DRA size. In the simulations we have been using the optimized N T /N U value reported in Table 4. We observe that the throughput dependency on N FR is rapidly decaying reaching an almost asymptotic performance for N FR > 20. This is a good result as it shows that the RRM performance is not very sensitive to the number of frames, hence the latency can be easily contained. When not stated otherwise, for the following simulations we adopted the value N FR = 30 which provides almost asymptotic throughput performance with an acceptable superframe duration. For this number of frames the throughput variation across the different frames for traffic case 2f was assessed to be less than 1 %.

4) Impact of the Number of Users/Slice Served
As explained in Sec. III, the proposed H-RRM algorithms are dividing the users to be served in the different frames composing the superframe, selecting for each frame a subset of users which can be mapped to beams maintaining a certain level of interference isolation from/to the other active beams. The number of active beams/frame allocated by the H-RRM algorithms is a random variable with an optimum average value N U which is dependent on the traffic distribution parameters. To exemplify the randomness of the number of users N U allocated on each frame by the H-RRM algorithm we report in Fig. 18 the simulation result for five distinct realization of traffic case 2f. Each trial histogram has a dif-  ferent colour to ease the reading of the results. It is apparent that the H-RRM is allocating a slightly variable number of users/frame (±5%) around the average value N U .
The most important parameters affecting the N U value are the number of antenna array elements N T and the traffic over overall potential satellite coverage area S T /S A ratio. As reported in Sec. IV-C5, we experimentally found that once the optimum N T /N U ratio has been derived for configuration 'a', its value can be extrapolated to another traffic case 'x' simply scaling it up with the ratio traffic areas ratio S x T /S a T as Fig. 19 shows. This is the approach we used to derive the optimum N T /N U ratio column in the following section Tables 4-5.

5) Impact of the non-uniform traffic cases
Tables 4-5 summarize the Shannon throughput simulation results for the various traffic configurations described in Sec. II-C. For each traffic configuration several simulations have been run to find the value of N T /N U which optimize the throughput. It can be remarked that the [N T /N U ] opt predicted value following the approach described in Sec. IV-C4 is working pretty well at least for the cases reported in Table  4.
Concerning the throughput performance for the various M-MIMO schemes with H-RRM we remark that for case 1 (full visible Earth uniform traffic distribution or PPP Distribution) the proposed H-RRM QVD approach is 13.8 % lower than the case 0 corresponding to the optimum RRM based on the ideal, yet unfeasible, minimum distance guaranteed PDP distribution.
It is remarkable to observe that QVD throughput is generally better than MB but also MF in several cases. The exceptions are represented by cases 2a, 2b (MF and MB) and 3a (MF only). Based on these findings the use of the QVD RRM algorithm is recommended with performance advantage as high as 9 %. The best performance are obtained for the case of even spots traffic probability reported in Table  5. However, the throughput difference between uneven and even spot traffic probabilities is limited to about 10 % which is considered negligible. This finding demonstrates the H-RRM robustness to the variation of traffic distribution. This is also confirmed by the limited difference in throughput between cases 2a vs 2b and 3a vs 3b. In the "b" cases σ XX is four times larger than "a" cases, causing the traffic PDF to appear much flatter (see Fig. 3-5-4-6). Nonetheless, the throughput reduction between 'a" and "b" cases is less than 7 %.
We also included the results for ideal MMSE M-MIMO which represents a performance upper bound. The MMSE throughput performance advantage is spanning from 7 to 12 % depending on the specific case. It is clear that, even disregarding the MMSE implementation complexity, in practice its performance will be well below to the reported values due to practical implementation aspects extensively discussed in [24].
It should be recalled that all the simulations have been performed with constant RF power, hence when the traffic area S T reduces the on-ground received power flux density increases accordingly. This means that the system is becoming more interference limited, hence, the MMSE potential advantage is slightly increasing as it is confirmed by the Sec. IV-C7 simulation results. Fig. 20 shows the simulated Shannon QVD throughput versus the ratio traffic area/coverage area for cases 2 and 3. The almost proportional throughput dependency to the ratio S T /S A is apparent, although the slope is more pronounced for case 3 than for case 2. It is evident that there is limited throughput difference between the cases of even and not even spot traffic probabilities which are quite different in terms of traffic density (see Fig. 3 and 4 for uneven spot probabilities and Fig. 5 and 6     tion results for the various traffic configurations described in Sect. II-C. As expected, the throughput is lower than for the Shannon case, but the results are in line with the previous findings.

6) Impact of Heuristic RRM Algorithm
So far the performance of the proposed pragmatic MB M-MIMO with H-RRM and QVD beam assignment scheme have been compared to an ideal MMSE precoding exploiting the same RRM scheme. To fully appreciate the importance of the proposed H-RRM solution we have been performing simulations using a simple RRM scheme which randomly assign the user to TDM frames with equal number of active users per frame. The beamforming/precoding matrix is then evaluated according to the population of users in the assigned colour. 9. In this case the H-RRM+QVD vs S-RRM+MMSE gain is up to 53 % for case 1 and 3d instead of a loss of -8 or -11 % when comparing H-RRM+QVD to H-RRM+MMSE.

7) Impact of RF power
As mentioned before, when reducing the traffic area served for constant payload RF power P c T the on-ground power flux density approximately increases proportionally to the ratio S A /S T . This means that the co-channel interference is also increasing. It is therefore interesting to assess the impact of a P c T reduction when the ratio S T /S A < 1. This has been investigated for the traffic case 2e where compared to case 2a there is a 70 % reduction in the traffic area. Tables 10 and 11 summarize the Shannon and DVB-S2X throughput simulation results. It is noted that reducing the RF power from 2.0 to 1.5 kW has a modest impact on the throughput (about 7 %). For P c T = 1.0 kW the throughput loss rises to about 17 %. Fig. 21 shows the simulation results for QVD and MMSE in a graphical format. It is apparent that the MMSE to QVD gap is reducing with P c T for the reasons explained before.  causing an increased co-channel interference mitigated by the MMSE. Anyway, this pragmatic M-MIMO QVD performance loss will be turned in a QVD gain when real MMSE implementation aspects will be considered. In addition the MB+QVD complexity is much lower compared to MMSE solution. Fig. 22 shows the simulated Shannon QVD throughput versus the ratio traffic area/coverage area for cases 2 and 3. The almost proportional throughput dependency in the ratio S T /S A is apparent, although the slope is more pronounced for case 3 than for case 2. Tables 14 and 15 summarize the Shannon and DVB-S2X throughput simulation results for the throughput dependency as a function of the payload RF power for the traffic case 2e. It is noted that reducing the RF power from 2.0 to 1.5 kW has a modest impact on the throughput (about 7%). For P c T = 1.0 kW the throughput loss rises to about 17%. Fig. 23 shows the simulation results for QVD and MMSE in a graphical format. It is apparent that the MMSE to QVD gap is reducing with P c T for the reasons explained before.    A will contain exactly n points is given by the Poisson's exponential function

V. SUMMARY AND CONCLUSIONS
We are interested in the probability distributionĜ(r) (known as Diggle's G-function [37]) of the distance D from a random point (event) to its nearest neighbor point (event) in the same point processĜ Exploiting the rule of subtraction of probabilities 9 we can interpret P (D < r) as the probability that a circle of radius r, centred on an arbitrary point contains at least one point, and P {D > r} as the probability that a circle of radius r, centred on an arbitrary point contains no point, which can be obtained substituting for n = 0 and A = πr 2 in the Poisson distribution P N (πr 2 ) = 0 = e −πρr 2 .
(64) 9 The probability that event will occur is equal to 1 minus the probability that event will not occur.
We thus obtain for the probability distribution of PPP which differentiated gives a Rayleigh probability density functionĝ The mean nearest neighbor distance of the PPP,R PPP can be evaluated as expectation of the Rayleigh probability density function The Poisson Disk Process, also known as Poisson Hard-Core Process (PHCP) [35], is a point process with a uniform distribution of occurrence of the points in the region of interest with the additional constraint that each point of the pattern must be separated from the others by a minimum distance 2r min , where the parameter r min is called the Poisson disk radius [39]. It is worth noting that in [24] it has been empirically shown that a user distribution according to PDP characteristics provides optimal throughput performance.
At the extreme end of the family of point patterns we have periodic lattices which are fully deterministic [38] . Given a point density ρ, the maximum mean nearest neighbour distance is obtained for hexagonal lattices

B. EMPIRICAL ANALYSIS OF POINT PATTERNS
The algorithms described in these technical note manipulate an input point pattern in a single (as in the case of VDI and QVD) or in multiple (as in the case of RRM) output point patterns. In all cases the input point pattern represents the geometrical distribution of users which is simulated as the realization of a set of random variables based on a continuous bi-dimensional probability density function. In first order the performance of the algorithms can be assessed with an analysis of the nearest-neighbor distance statistics of the output point patterns. For this reason we need to discuss/adapt some of the techniques needed for point pattern analysis.
For each point of the point pattern under analysis we can evaluate the nearest-neighbor distance and build the empirical distribution function, a staircase function with steps of height 1/N (where N is the number of points) for each of the ordered nearest-neighbor distances.
In plotting the empirical distribution function of the nearest-neighbor distance it must be considered that the distance is an absolute geometric measure and to compare different distributions we need to define a re-normalization of the x-axis (the y-axis is intrinsically normalized to 0-1 by the definition of distribution function).
Considering the centrality of the Poisson point process for point pattern analysis (including hypothesis testing of the complete spatial randomness) we can introduce a re-scaling of the distance to the PPP mean nearest-neighbor distance for the given point density. The re-scaling procedure follows: • Given A, the area of the Region of Interest, and N , the size of the empirical point pattern population, define the empirical density • Re-scale to distances, r, to the PPP mean nearestneighbor distance for the empirical density ρ r = r R P P P (ρ) = 2r The re-normalization allows to compare point patterns with different sizes of the number of empirical points. Examples of point patterns and relevant empirical distribution function of the nearest-neighbor distance are reported in the following  ity density function (PDF) given the empirical observations [40]. It is based on the optimization of the PDF parameters such that the likelihood of making the observations is maximised for the the optimal parameters. MATLAB statistical toolbox has a generic MLE routine, mle( ), which has been used for the estimations of the of the RRM Diggle's Gfunction.  Concerning the functional form of the PDF for the best fitting, a Weibull three-parameters PDF [45] has been implemented (being not available in MATLAB).
(71) Its adoption is empirically justified on the base of the following elements: • Considering that the Poisson Point Processes has an exact theoretical Rayleigh distribution, which is a oneparameter (i.e. α -scale) distribution, it is natural to consider the Weibull distribution which is a two-parameter generalization (i.e. α -scale and β -shape) of the Rayleigh distribution (for β = 2 the Weibull coincides with Rayleigh distribution) . • The additional shape parameter of the Weibull distribution allows to change the growth rate (for an β > 2) of the PDF allowing to adapt it to different empirical slopes.  • The third parameter (r 0 -location) added to the standard Weibull distribution allows to control both the non-zero attack point of the distribution, and to adapt a given slope to the mean. The qualitative goodness of the estimation for a typical RRM point pattern can be appreciated in Figure 26 (compare the black empirical PDF with red Weibull 3P best fitted with MLE techniques).