Base Station Placement Algorithm for Large-Scale LTE Heterogeneous Networks

Data traffic demands in cellular networks today are increasing at an exponential rate, giving rise to the development of heterogeneous networks (HetNets), in which small cells complement traditional macro cells by extending coverage to indoor areas. However, the deployment of small cells as parts of HetNets creates a key challenge for operators’ careful network planning. In particular, massive and unplanned deployment of base stations can cause high interference, resulting in highly degrading network performance. Although different mathematical modeling and optimization methods have been used to approach various problems related to this issue, most traditional network planning models are ill-equipped to deal with HetNet-specific characteristics due to their focus on classical cellular network designs. Furthermore, increased wireless data demands have driven mobile operators to roll out large-scale networks of small long term evolution (LTE) cells. Therefore, in this paper, we aim to derive an optimum network planning algorithm for large-scale LTE HetNets. Recently, attempts have been made to apply evolutionary algorithms (EAs) to the field of radio network planning, since they are characterized as global optimization methods. Yet, EA performance often deteriorates rapidly with the growth of search space dimensionality. To overcome this limitation when designing optimum network deployments for large-scale LTE HetNets, we attempt to decompose the problem and tackle its subcomponents individually. Particularly noting that some HetNet cells have strong correlations due to inter-cell interference, we propose a correlation grouping approach in which cells are grouped together according to their mutual interference. Both the simulation and analytical results indicate that the proposed solution outperforms the random-grouping based EA as well as an EA that detects interacting variables by monitoring the changes in the objective function algorithm in terms of system throughput performance.


Introduction
The last two decades have witnessed a boom in the use of cellular communication technologies. Billions of people are now requesting high-quality mobile wireless services with end-user data the minimum set of BSs is selected from a fixed set of candidate sites that satisfy quality and outage constraints. The placement of BSs is then determined in a subset of the deployment area according to private property limitations or electromagnetic radiation constraints. However, while the site-selection and site-placement algorithms provide a more locally optimal solution with a lower number of BS, a larger number of BSs requires more computation time.
Because evolutionary algorithms (EAs) can be characterized as global optimization methods, they have been utilized successfully in a variety of complicated real-world applications [13]. Several methods have also been developed spontaneously in the field of radio network planning, mainly based on EAs [14][15][16]. Weicker et al. [14] proposed a steady-state EA with Pareto tournaments (stEAPT), which considers frequency assignment and channel interference for BS placement. Most recently, addressing the observation that multiple objectives (MO) must be taken into account when solving the wireless heterogeneous transmitter placement problem, Ting et al. [15] proposed to integrate a novel variable-length representation and a new crossover approach into their non-dominated sorting genetic algorithm II (NSGA II) [16], which is known for its effectiveness in dealing with MO problems. However, one crucial difficulty in employing EAs is the huge time consumption resulting from the high complexity of performance analyses for fitness evaluation and the large number of evaluations needed in evolutionary optimization techniques. Accordingly, the performance of EAs often deteriorates rapidly with the growth of search space dimensionality.
To overcome the problems mentioned above when designing optimum network deployments for large-scale LTE HetNets, we attempt to decompose the high-dimensionality problem and tackle its subcomponents individually. We propose a grouping method to divide the candidate solutions (individuals) in the populations into groups. Noting that some HetNet cells have strong correlations due to inter-cell interference, we propose to use a correlation grouping approach instead of grouping the individuals randomly with the aim of rapidly converging to optimal solutions. In this approach, variables with strong correlations (i.e., interfering cells) form a group when finding the optimal deployment of heterogeneous cells in the HetNet. In addition, we modify the variable-length genetic algorithm presented in [15] to be applied to the divided groups.
The rest of this paper is organized as follows: In Section 2, we include a short review of how EAs and their variants have been applied for solving the problems in various applications including the BS placement problem. Section 3 presents the mathematical formulation of the BS deployment problem optimization in the LTE HetNet. Section 4 presents the proposed grouping method and the strategy for solving the optimization problem based on the variablelength genetic algorithm [15]. In Section 5, we provide an analytical model of the probability of two BSs placed in the same group interfering with each other, as well as the corresponding numerical results, and report simulation results for different user distributions. Finally, conclusions are drawn in Section 6.

Related Work
EAs are characterized as global optimization methods and are generally known to be robust optimizers that are well suited for objective functions that are discontinuous and have many non-smooth changes [17,18]. For this reason, they have been applied successfully to a variety of complicated real-world applications such as discovery of link communities in complex networks [13], financial and economic applications [19], aircraft conflict avoidance [20,21], demand side management in smart grids [22], etc.
The BS placement problem is to find the optimal positions of BSs, considering various controlled and uncontrolled factors of traffic density, capacity, interference, existing BSs, etc. Due to the combined effects of these factors, the problem cannot be solved in polynomial time; that is known to be NP-hard [4]. Some heuristic methods based on the evolutionary paradigm have also been developed for finding high-quality solutions to such BS placement problems [14,15]. In [14], the steady-state EA with Pareto tournaments (stEAPT) approach is introduced as a new MO technique that considers frequency assignment and channel interference for BS placement. This approach combines a steady-state scheme with a very efficient data structure leading to superior time complexity. Most recently, addressing the observation that MO must be taken into account when solving the wireless heterogeneous transmitter placement problem, Ting et al. [15] proposed to integrate a novel variable-length representation and a new crossover approach into their non-dominated sorting genetic algorithm II (NSGA II) [16], which is known for its effectiveness in dealing with MO problems. However, one crucial difficulty in employing EAs is the huge time consumption resulting from the high complexity of performance analyses for fitness evaluation and the large number of evaluations needed in evolutionary optimization techniques. Accordingly, the performance of EAs deteriorates rapidly with the growth of search space dimensionality. Apparently, this is also the case with the BS placement problem in LTE HetNets because explosive mobile data demands have driven mobile operators to deploy LTE small cells on a large scale.
Cooperative co-evolution (CC) has been introduced into EAs with the aim of solving increasingly large and complex optimization problems through a divide-and-conquer paradigm [23]. Nonetheless, existing CC algorithms did not take into account variable interdependencies for nonseparable problems in which tight interactions exist among different decision variables. To efficiently tackle nonseparable problems, some CC frameworks were proposed relying on random grouping that randomly allocates the variables to subcomponents in every co-evolutionary cycle [24][25][26], instead of using a static grouping. These algorithms do not provide a systematic procedure to group the interacting variables nor to detect their interdependencies, even though it was shown in [24] that with random grouping, the probability of placing two interacting variables in the same subcomponent for several cycles is reasonably high.
More recently, some algorithms have been proposed to identify and group interacting variables into common subcomponents in various real-world optimization problems [20,21,23]. In [20] and [21], CC frameworks with dynamic grouping strategies were proposed to guarantee safety in air traffic control. In the dynamic grouping strategy, a large number of aircraft are divided into several sub-groups based on their interdependence and the sub-groups are adjusted dynamically as new conflicts appear after each iteration. Omidvar et al. [23] proposed a decomposition method called differential grouping that is able to group the interacting variables with high accuracy, focusing on large-scale global optimization problems. Here, it should be noted that these algorithms focus on discovering the interdependencies between variables, whereas in the LTE HetNets, the effect of inter-cell interference should be designed as various levels of interdependence, since the interference is one of the most critical factors to be considered when deploying small cells. Specifically, the degree of interdependence between cells varies based on the amount of inter-cell interference, not the existence of interdependence. Thus, we cannot efficiently group the cells with strong interdependencies by simply applying the existing grouping methods to the BS placement problem.
In this paper, we propose to use a correlation grouping approach with the aim of rapidly converging to optimal solutions. In this approach, variables with strong correlations (i.e., interfering cells) form a group when finding the optimal deployment of heterogeneous cells in the HetNet. In addition, we modify the variable-length genetic algorithm presented in [15] to be applied to the divided groups.

BS Deployment Optimization Problem Formulation
In this section, we present a formulation of the BS deployment optimization problem in the LTE HetNet, defining the objective, variables, and constraints for the problem. In this study, we are given a domain, D, in the HetNet that must be covered by heterogeneous BSs with different transmission powers. For the given domain, we aim to find a BS deployment plan that maximizes user satisfaction in terms of the throughput provided per unit of traffic demand from the user.
We assume that there are M BSs in the domain D, and that each BS is located at the position (x b , y b ) 2 D (1 b M) with the transmission power p b , where M is a constant value determined by the network operator. Let p max be the maximum transmit power and 0 p b p max . It is also assumed that there are N users and that each user is located at the position (X u , Y u ) (1 u N) with a traffic demand, d u .
The signal-to-interference-plus-noise-ratio (SINR) at the u-th user is then given by where N 0 is the noise power and l b,u is the path loss between the BS, b, and the user, u, which is [27].
We denote with T u and U u the LTE downlink throughput and the user satisfaction for u-th user, respectively. Given that the u-th user is served by the b-th BS, we can express T u and U u as follows: and where R u is the number of resource blocks (RBs) assigned to the u-th user, and C = 180kHz is the bandwidth of an RB [27]. For a given ordered set of triples containing the position and the traffic demand of every user, the system satisfaction function F U is defined as for the ordered set of triples containing the position and transmit power of every BS, Let A be the set of all possible BS deployments in the domain D with a maximum transmission power of p max , and let H b be the set of all users served by the b-th BS. The problem of finding the optimal BS deployment in A is then formulated as follows: where R max is the maximum number of RBs that can be allocated, being set to 50 when using a 20 MHz system bandwidth in LTE frequency division duplex (FDD) downlink [28].

Cooperative Co-evolution with BS Grouping
In this section, we describe the proposed grouping method and the EA used to solve the BS placement optimization problem.

BS Grouping Algorithm
As addressed in Section 1, increased wireless data demands have driven mobile operators to roll out large-scale networks of LTE small cells. Accordingly, the problem domain presented in the previous section is considered quite large, leading to an increased complexity in solving the problem. We note that cooperative co-evolution has been proposed to solve large and complex problems through problem decomposition [29]. Based on this notion, we decompose the BS placement optimization problem into subproblems by dividing all the BSs into different groups.
The subset B j is obtained as follows: 1. G pivot points are randomly chosen in the domain D, each of which represents each group B j .
2. Each BS is included in a subset B j , in which the BS has the highest received signal strength among G disjoint subsets.

Proposed Evolutionary Algorithm
In this section, we present the details of our proposed EA for finding the optimal solution to the BS placement problem. Following the general concept of EAs, the algorithm starts with a population P composed of a set of individuals B (i) . The proposed EA is iterated until the solution converges to the optimal placement. Each EA consists of the following four steps: fitness evaluation, parent selection, crossover and mutation. 1) Fitness evaluation. At the beginning of each iteration, all the individual's fitnesses are evaluated by computing the objective function (i.e., fitness function) in Eq (5). Thus, the higher the system throughput, the higher the fitness value. Then, only the best 50% are retained for the next generation.
2) Parent Selection. For each of the groups generated by the BS grouping algorithm presented in Section 4.1, a pair of parents is selected for the next generation. Specifically, for each group in the individual B (i) , B ðiÞ j itself becomes a parent, while the other parent B ði 0 Þ j (i 0 6 ¼ i) is selected from the best 50% of individuals, excluding the parent already selected, with the probability If the typical parent selection strategies available in the literature [14,15] were applied to our EA, it would not be possible to decide which of its two parents a child belongs to. However, no such ambiguity exists in the parent selection process described above since only one child is generated for each individual.
3) Crossover. A child individual is produced from the two selected parents through crossover and mutation. Given the crossover rate P c , the hybrid crossover has three ways of producing a child [15]: 1) perform only uniform crossover with the probability P c × P c ; 2) perform both uniform crossover and one-point crossover with the probability P c × (1 − P c ); or 3) perform only one-point crossover with the probability (1 − P c ).
Regarding the representation of the position and transmission power of BS, B (i) , (x i , y i , p i ) consists of 24 bits, where the leftmost 16 bits (2 × 8 bits) and the rightmost 8 bits indicate the position and the transmission power, respectively. Then, we apply the uniform crossover method equiprobably to choose which of the two parents the child will inherit a bit from.
We now propose to modify the one-point crossover between two parents, B ðiÞ j and B ði 0 Þ j . We denote by j B ðiÞ j j the number of tuples consisting of B ðiÞ j , that is the number of BSs belonging to the set, B ðiÞ j . Then, the proposed one-point crossover proceeds as follows: Step 1). A random point is chosen in the range of ½0; minðj B ðiÞ j j; j B ði 0 Þ j jÞ. Each parent is divided into two parts at this point.
Step 2). A child is produced by combining the first part of B Under the assumption made regarding M, the total number of BSs, M, in the domain, D, is constant. However, while the length of a child is variable in [15], leading to a value different from M, the one-point crossover proposed here maintains the length of the individual to be M over subsequent generations. Examples of the one-point crossover presented in [15] and the proposed method are depicted in Fig 1. It can be seen in Fig 1 that the two children have different lengths (7 and 4 tuples) from those of Parent 1 and Parent 2 (5 and 6 tuples, respectively) in the one-point crossover [15], while the proposed EA keeps the length of the child the same as that of its parent. 4) Mutation. After the crossover, a bit-flip mutation with the mutation rate P m is performed, giving each of the child's bits a chance to flip.

Results
As explained in the Section 1, the proposed correlation grouping approach, which enables interfering cells to be placed together in one group, is a key contribution, since some cells have a strong correlation due to inter-cell interference in the HetNet. Thus, in this section, the probability of placing two interacting variables into a single group is first analyzed numerically for our proposed approach and a random grouping scheme [24][25][26]. We then provide evidence, based on simulation results, in support of throughput performance improvement of the proposed scheme over the random grouping scheme.

Probability of Interacting Variables Belonging to the Same Group
We denote by P g the probability of two BSs being placed together in a single group at least N k times during N c cycles, where a cycle consists of one complete evolution of all groups. We simply refer to this probability as the grouping probability in this paper. In the random grouping scheme, the grouping probability P g is derived as follows [24][25][26]: Note that the SINR is subjected to considerable attenuation over distance between the transmitter and the receiver as seen in Eq (1). Thus, we derive the grouping probability in our proposed EA based on the distance between BSs. We assume that BSs are distributed over the domain D according to a Poisson point process (PPP) with mean density, l ¼ M AðDÞ , where A(D) denotes the area of the domain D. It is also assumed that all BSs use the same transmission power and have circular coverage. Then, we compute the probability that two BSs, b i and b i 0 are placed in the same group for the following two cases: Case 1) Either of b i or b i 0 is chosen as the pivot point; Case 2) Neither b i nor b i 0 is chosen as the pivot point.
Case 1). Let d i,i 0 be the distance between b i and b i 0 . If b i is chosen as one of G pivot points, then the other BS, b i 0 must belong to the same group as b i . Thus, none of the G pivot points minus the chosen pivot b i 0 can be located within a circle of radius d i,i 0 centered at the BS b i since the interferer's transmission power is known to depend on the distance to the BS, b i 0 , as can be seen in Eq (2).
Let f M,G (n) be the probability of selecting n BSs among M BSs minus G pivot points. We then calculate the probability f M,G (n) as follows: Using Eq (11), the probability of selecting the BS b i 0 given that b i is the pivot point, P(iji 0 ) is (x j 1 , y j 1 , p j 1 ) ... (x j n j , y j n j , p j n j ) Evolutionaly algorithm in [13] 3 BSs 4 BSs given by where p i 0(d i,i 0, n) is the probability that a circle of radius d i,i 0 centered around b i 0 contains exactly n points, which is given by following the PPP model. Then, in Case 1, we can derive the probability that two certain BSs, b i and b i 0 , belong to a single group, P ð1Þ c , as follows: where the first term, (G(M−G))/(M(M − 1)) indicates the probability of Case 1 occurring. Case 2). Consider an arbitrary BS and suppose that this BS is the n-th nearest one to the BS b i (let's say b i 00). We are now interested in obtaining the probability that the BS b i belongs to the group represented by b i 00, which is denoted by P(iji 00 ). Note that the probability P(iji 00 ) is equivalent to the probability that the BS b i 00 becomes the nearest one of b i . Given the distance d i,i 00 from the BS, b i to the n-th nearest BS, it follows that the probability P(iji 00 ) is given by , using the PPP model.
We next derive the probability that the BS b i 0 also selects b i 00 as its pivot point P(i 0 ji 00 ), conditioned on the fact that the BS b i belongs to the group represented by b i 00. Let A i (r) and A i 0 (r) denote circular areas with a radius of r centered around the BSs b i and b i 0 , respectively. Then, we can express the probability P(i 0 ji 00 ) as follows: wherep i 0 ðD i 0 ;i 00 ; nÞ is the probability that a subdomain D i 0 ;i 00 ¼ D \ fA i 0 ðd i 0 ;i 00 Þ À A i ðd i;i 00 Þg contains exactly n points. Given the occurrence of Case 2, from Eqs (14) and (15), we get the probability that both b i and b i 0 belong to the group represented by a pivot point b i 00, P ð2Þ c as follows: where the first term indicates the probability of Case 2 occurring.
Finally, from Eqs (13) and (16), we can derive the grouping probability in our proposed scheme as follows:

Numerical Results
The results presented in [24] confirmed that their random-grouping approach performs better than EAs without grouping when tackling large optimization problems. Thus, in order to assess the accuracy of our analysis in terms of grouping probability, we carried out simulation tests for both the random and proposed grouping schemes, and compared the predictions of the analytical model from the previous section with these simulations. It can be seen in Fig 2(a) that both the numerical and the simulated grouping probabilities decrease as d i,i 0 increases in the proposed grouping scheme. This is because inter-cell interference decreases with an increase in the distance between the two BSs. On the other hand, as seen in Fig 2(b), all the curves for the random grouping are flat in the entire range of d i,i 0 since the random grouping scheme does not consider inter-cell distance.
When producing the results for N c = 30, N k is set to 15, which indicates the probability that two BSs are assigned to one group for at least 15 cycles. The grouping probability in the proposed grouping decreases quickly when d i,i 0 rises above 400 m and 300 m for G = 5 and 9, respectively, and the decrease for G = 9 is more rapid than that for G = 5. Even for G = 2, the grouping probability is flat regardless of the change in inter-cell distance. The reason for this phenomenon is that the smaller the number of groups, the lower the probability that two BSs belong to the same group. We also observe from Fig 3(b) that the grouping probability is the same regardless of the value of d i,i 0 when random grouping is employed. As can be seen in Figs 2 and 3, the predicted and the observed results correspond closely, with only a slight difference between the two.

Simulation Environment
To evaluate the advantages of the proposed EA by means of simulation, we distributed users in a square area of 2 km × 2 km. To simulate an urban area with high user density [30], the user density was set to 10 users/km and 15 users/km, which are equivalent to 400 and 900 users, in the entire simulation area, respectively. Note that the user distribution affects the amount of traffic demand in an area, which should be considered when placing the BSs as presented in Section 3. To evaluate the effectiveness of the proposed EA in different user distributions, we used three user distribution models in the simulations: a uniform distribution, a Gaussian distribution, and a four-Gaussian hotspot distribution [4]. The three different user distributions are illustrated in Fig 4. In the uniform distribution, users are distributed uniformly over the entire area. The Gaussian distribution models a hot spot with user density at a maximum located at the center, and gradually decreasing toward the boundary. The four-Gaussian hotspot distribution models four hot spots with very densely located users. We omit the illustrations of the three distributions for the case of 15 users/km because, except for the density, they are the same as those shown in Fig 4. All users were assumed to have the same traffic demand of 1 Mbps. With regard to the number of BSs, M was set to N/10, which is 10 users per BS. In the simulation, the value of the mutation rate is configured using P m = 1/substring_length, where the substring length is 24 and the crossover rate is set to 0.9 [15]. As suggested in [23] and [31], the population size is set to 50. We set the maximum transmission power, p max , to 46 dBm [27,30]. The simulation parameters are summarized in Table 1.

Simulation Results
As mentioned in Section 3, we started initially with M randomly located BSs, which will be referred to simply as 'M-Random' and then ran the proposed EA presented in Section 4 to find the optimal locations of the BSs. The simulation was performed in two scenarios, without and with macro BSs being installed initially. We first compare the performance of the proposed EA with that of M-Random to demonstrate that large gains cannot be achieved merely by installing more BSs.    It is shown in Fig 5 that in M-Random, the slope of the increase tapers off even though the system throughput increases as M increases. This is because deploying too many BSs randomly can cause severe interference; that is to say, it is more likely for a user to suffer severe interference from a neighbor cell covering the same area as the serving cell of the user. Notably, the BS deployment that is obtained from the proposed EA starting with M = 40 (M = 90) achieves a system throughput almost as high as in M-Random with M = 80 (M = 120). That is, the proposed EA is able to achieve the same throughput by deploying a much smaller number of BSs than the random deployment approach. We can also see that the proposed EA shows a higher system throughput for the two Gaussian distributions than that for the uniform distribution, whereas the contrary phenomenon is observed in M-Random because the proposed EA finds the BS deployment that maximizes system throughput, as can be seen in Eq (7).

System Throughput Improvement
We now present the simulation results for the proposed EA, compared to the random-grouping based EA (RG-EA) [24][25][26] and an EA that detects interacting variables by monitoring the changes in the objective function and groups them as in [23]. We name the latter algorithm as EA with grouping based on interaction-detection (IDG-EA). When simulating IDG-EA, the threshold to identify the interacting variables is set to 10 −3 , as recommended in [23].
Because we aim to maximize the system throughput per user (i.e., F U ) in this study, the system throughput improvement (denoted by ΔF U ) is used as the main metric to evaluate the performance of the three algorithms. More specifically, ΔF U indicates the amount of throughput improvement (in percentage) by the three algorithms over the initial random deployment of M BSs. Tables 2-4 show the average, standard deviation, minimum and maximum of ΔF U in RG-EA, IDG-EA, and the proposed EA for each of the three user distributions, without any macro BS being initially installed. Observe in Tables 2-4 that the proposed EA improves the performance of RG-EA by up to 18.97% and 20.68% for the two user densities of 10 and 15 users/km, respectively. Notably, when the uniform distribution is used and the user density is 10 users/km, the minimum value of ΔF U in the proposed EA is even greater than that of RG-EA. We also observe that the statistics of ΔF U in the proposed EA are higher than those in RG-EA for the other two user distributions. Specifically, ΔF U of the proposed EA is up to 16.64% and 12.33% (17.08% and 12.52%) higher than that of RG-EA for the Gaussian distribution (the four-Gaussian hotspot distribution) when there are 10 and 15 users/km, respectively. Tables 5-7 present the statistics of ΔF U in the case that BSs are already installed in the simulation of an urban area in which some BSs are already installed such that inter-site distance (ISD) is 750 m [32]. For all three user distributions, the proposed EA shows a larger ΔF U than that of RG-EA for user densities of both 10 and 15 users/km. Specifically, the average values of ΔF U in the proposed EA are up to 11.78%, 16.28%, and 8.89% (10.89%, 11.57%, and 7.50%) higher than those of RG-EA for the uniform, the Gaussian, and the four-Gaussian hotspot user distributions, respectively, when the user density is 10 users/km (15 users/km).
We see from the six tables that the proposed EA also achieves higher system throughput than IDG-EA. Specifically, when no macro BS is initially installed, the average values of ΔF U in the proposed EA are up to 54.76%, 73.41%, and 104.47% (71.04%, 54.06%, and 94.42%) higher than IDG-EA for the uniform, Gaussian, and four hotspot user distribution, respectively, when the user density is 10 users/km (15 users/km). In the case that some BSs are installed, the average values of ΔF U in the proposed EA are up to 37.34%, 69.02%, and 32.38% (43.18%, 89.84%, and 40.97%) higher than IDG-EA for the uniform, Gaussian, and four hotspot user distribution, respectively, when the user density is 10 users/km (15 users/km). We also observe that RG-EA outperforms IDG-EA as well. This is because in IDG-EA, more separable variables are classified as nonseparable ones. Interestingly, we also find in the tables that maximum system throughput improvements are achieved for M = 20 or 30 when there are 10 users/km, while M = 30 or 40 when there are 15 users/km, and that after the point at which the maximum occurs, the throughput improvement starts to decrease. This is because installing more BSs causes harmful inter-cell interference, leading to a reduction in the system throughput improvement. Finally, we can say that the proposed EA outperforms both RG-EA and IDG-EA, not only for the uniform user distribution, but also for the other two practical user distributions.

Conclusion
In this paper, we have presented a correlation grouping approach for the application of an EA with the aim of designing an optimum network planning algorithm for large-scale LTE Het-Nets that results in the rapid convergence to optimal solutions. Noting that some HetNet cells have a strong correlation due to inter-cell interference, the correlation grouping approach makes variables with strong correlations (i.e., interfering cells) form groups, instead of grouping the individuals randomly when finding the optimal deployment of heterogeneous cells. We have also modified the variable-length genetic algorithm presented in [15] to be applied to the divided groups. To evaluate the performance of the proposed algorithm, we have analyzed the grouping probabilities and conducted simulations. Both numerical and simulation results confirm that the proposed algorithm outperforms both RG-EA and IDG-EA in terms of system throughput, not only for the uniform user distribution, but also for the practical Gaussian and four-Gaussian hotspot user distribution models.