Two-Step Dynamic Cell Optimization Algorithm for HAPS Mobile Communications

High altitude platform stations are attracting much attention as novel mobile communication platforms for ultra-wide coverage areas and disaster-resilient networks. Multi-cell configurations can increase communication capacity in a wide area. Conventional work optimizes the cell configuration for a multi-cell configuration based on a genetic algorithm (GA). This method identifies an optimal cell configuration in terms of spectral efficiency depending on the number of cells under a uniform user distribution. However, user distributions differ depending on location; thus, cell configuration optimization is also required for non-uniform user distributions. The problem is an increased number of parameters because each cell requires different antenna parameters for non-uniform user distribution compared to a uniform user distribution. This increases the difficulty of optimization, even when using a GA in some cases. Therefore, we propose a GA-based two-step cell optimization algorithm that comprises both search space reduction and antenna parameter optimization steps to address this problem. The proposed method employs the concept of co-evolution, i.e., a divide-and-conquer method. The proposed method divides multiple cells into several groups to reduce the number of optimized parameters and optimizes each group in order. In addition, the search space is reduced based on a marginal histogram obtained via optimization with a large step size. Simulation results demonstrate that this technique can reduce the number of combinations compared to the case without sub-area division. In addition, there are cases where the search space reduction further reduces the number of combinations without degrading the sum of the square root of throughput.


I. INTRODUCTION
High altitude platform stations (HAPS) offer great potential for innovation in the mobile communication services field [1], [2]. HAPS mobile communications are expected to provide mobile communication services directly to smartphones on the ground, which are commonly used in terrestrial mobile communication networks with fourth-generation long-term evolution (4G LTE) and fifth-generation new radio (5G NR). HAPS mobile communications employ unmanned aerial vehicles (UAV), e.g., balloons and airships that fly in the stratosphere at an altitude of 20 km, where low wind speeds are expected. The altitude of HAPS is significantly less than that of communication satellites, e.g., 36,000 km The associate editor coordinating the review of this manuscript and approving it for publication was Mohammad Tariqul Islam .
for geostationary earth orbit (GEO) satellites, 2,000-36,000 km for medium earth orbit (MEO) satellites, and less than 2000 km for low earth orbit (LEO) satellites. The round trip time (RTT) from user equipment (UE) to a satellite is 240-280 ms for GEO satellites, 54-86 ms for MEO satellites, and 6-30 ms for LEO satellites [3]. In contrast, the RTT from UE to a HAPS is 360 µs and 680 µs at 50 and 100 km from the center, respectively. Due to its low latency characteristic, a HAPS is a promising platform for 5G to enable URLLC (ultra-reliable and low latency communication) where 3GPP defines a target packet failure rate of 10 −5 within 1-ms overthe-air latency [4].
The coverage area of a HAPS ranges from several dozen kilometers to 100 km in radius, which necessitates a multi-cell configuration with single-cell frequency reuse to increase communication capacity [5], [6]. By dividing a VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ coverage area into multiple cells, each cell can have a narrow beamwidth, which improves the antenna gain and link budget, thereby extending the coverage area. However, the optimal cell configuration in terms of spectral efficiency (SE) for a multi-cell configuration is unclear and is dependent on the number of cells. To clarify the optimal cell configuration, we previously proposed a cell configuration optimization method based on a genetic algorithm (GA) [6]. This method optimizes antenna parameters to maximize SE under the assumption that users are distributed uniformly in the coverage area. However, in practice, user distributions are not necessarily uniform depending on location. In addition, sudden increases in traffic demand can occur in specific areas, e.g., during a disaster. Thus, a cell optimization method that considers various user distributions is required. Previous studies have focused on multi-beam optimization and beamforming in satellite communications [7]- [14] and control beam position considering user location [7], [8]. For example, a previous study [7] proposed a dynamic coverage adjustment scheme for satellite mobile communications according to the spatial distribution of user traffic. However, this scheme only considers dynamic beam coverage switching between two kinds of coverage patterns. In another study [8], beams were optimized to attempt to cover all users with the smallest number of cells. In addition, the studies in [9], [10] presented beamforming schemes to maximize the sum-rate for satellite-terrestrial integrated networks and satellite and aerial-integrated networks, respectively. For terrestrial networks, several studies have investigated coverage and capacity optimization (CCO) of self-organizing networks. For example, a previous study [15] proposed a self-tuning sectorization method using deep reinforcement learning to optimize broadcast beams autonomously and dynamically based on user distribution. This method considers UE's measurement results as input and optimizes beamwidth and direction to maximize the number of connected users. However, these methods [8], [15] do not consider potential users that appear suddenly. In addition, other studies addressed CCO using reinforcement learning [16], [17]. However, it is difficult to learn every situation in advance, thereby making the output from reinforcement learning is non-optimal if learning is inadequate. Thus, a fast optimization method that can manage variable situations and user distributions is required.
GAs are candidate solutions for fast optimization because they are easy to parallelize. Although each cell can have the same antenna parameters among cells in a uniform user distribution, each cell requires different antenna parameters to realize optimization in non-uniform user distributions. As a result, the parameter combinations search space expands exponentially as the number of cells increases, thereby making it difficult to optimize parameters in some cases, even when using a GA.
Techniques to solve complicated problems with many parameters can be classified into problem decomposition and search space reduction techniques. In terms of problem decomposition, a complicated problem can be solved by a divide-and-conquer method, e.g., the co-evolutionary algorithm [18], which was inspired by mutualism, i.e., an ecological interaction between two or more species where each species has a net benefit. In a co-evolutionary algorithm, a complicated problem is decomposed into multiple subproblems, and each subproblem is solved individually.
To optimize the cell configuration for non-uniform user distributions, we propose a GA-based two-step cell optimization algorithm comprising a search space reduction and antenna parameter optimization steps. Based on the co-evolution concept, the proposed method divides multiple cells into several groups to reduce the number of parameters to optimize at a time and optimize each group in order. The proposed method optimizes antenna parameters, e.g., beamwidth and direction in a group while maintaining certain communication quality for the entire area. To further reduce the combinations of parameters, in the search space reduction step, optimization with a large step size is performed to identify promising regions from the original large search spaces. Here, the search space is reduced based on a marginal histogram of each parameter, and the reduced search range is explored with the desired step size in the optimization step. The proposed method is beneficial, for example, when many HAPSs are deployed to provide nationwide coverage. It is because configuring many HAPSs manually according to non-uniform user distributions requires lots of effort. In addition, tracking changes in the user distributions over time can be realized by executing this algorithm regularly.
We present simulation results to validate the proposed method. The results demonstrate that sub-area division is effective in terms of reducing the number of performed combinations and improving the value of an objective function, which is defined as the sum of the square root of throughput for all users. In addition, we demonstrate that there are cases when the number of performed combinations can be reduced without degrading the objective function value by performing the search space reduction step.
The remainder of this paper is organized as follows. Section II describes the system model. The proposed two-step dynamic cell optimization algorithm is presented in Section III. The genetic algorithm and problem formulation are presented in Section IV, and simulation results are presented and discussed in Section V. Finally, the paper is concluded in Section VI.

II. SYSTEM MODEL
In this section, we describe the system model of HAPS mobile communications, as shown in Fig. 1. A HAPS provides mobile communication services in a circular area with a radius of up to 100 km at most. In this paper, we consider the deployment of a single HAPS. The primary system parameters are listed in Table 1. The frequency is 2 GHz, which is identified as a HAPS mobile communication service band in the International Telecommunication Union Radiocommunication Sector (ITU-R). Note that the propagation model assumes free space path loss due to a good line of  sight environment. Following the literature, the antenna gain for the UE is assumed to be −3 dBi [20]. The noise figure is assumed to be 5 dB for the UE and 3 dB for the base station (BS), which corresponds to realistic device and BS implementations. The antenna gain depends on the beamwidth; thus, the downlink (DL) antenna gain at the transmitter and uplink (UL) antenna gain at the receiver are denoted ''Design parameters'' in Table 1. We set the bandwidth to 18 MHz, which is the largest bandwidth available in an LTE system. The UL bandwidth is 360 kHz, which corresponds to two resource blocks (RBs) and is required to achieve a UL connection.  Fig. 2(a) is a two-layer configuration comprising inner and outer layers, and the 21-cell configuration shown in Fig. 2(b) is a three-layer configuration. The parameters required to form a cell are vertical beamwidth θ 3dB , horizontal beamwidth φ 3dB , vertical direction angle θ str , and horizontal direction angle φ str for each cell, as shown in Fig. 3. Although each cell has four parameters, the parameters except horizontal direction in the same layer can have the same value in a uniform user distribution scenario. For a uniform distribution scenario, the horizontal direction angle of each cell can be spaced equally, and it is sufficient to designate an offset to rotate all cells in a layer. Thus, the total number of parameters in a uniform user distribution scenario is four multiplied by the number of layers, e.g., eight for the nine-cell configuration and 12 for 21-cell configuration. In this paper, we consider the cell configurations with N = 9 and N = 21.

B. ANTENNA PATTERN
The antenna pattern considered in this study is based on [6], which is modified from ITU-R recommendations [19]. The antenna gains for angle can be expressed as where b is the half-power beamwidth divided by 2, 1 = b √ −L N /3, 2 = 3.745 b , X = L N + A log 10 ( 2 ), and 3 = 10 (X −L F )/A . Based on a real antenna pattern that can be realized for a planar patch array antenna, we modify the near-in-side-lobe level L N , the far side-lobe level L F , and A to −20 dB, −30 dB, and 20, respectively. The antenna gain can be expressed as follows: where G v and G h are the antenna gains for the vertical and horizontal planes, respectively, as calculated in (1). Note that is omitted for the expression for G, G v , and G h . The maximum antenna gain G p is expressed as follows: which is based on an antenna with a 6 dBi gain and beamwidth of 80 • for both the vertical and horizontal planes.

C. COVERAGE CONDITIONS
To realize a certain communication quality in the entire coverage area, we define coverage conditions for the signal-tonoise ratio (SNR) and signal-to-interference-plus-noise ratio VOLUME 10, 2022 (SINR). Specifically, 99% of all users assuming uniform distribution in the area should satisfy the following constraints [6].
• DL SNR above 8 dB • UL SNR above 11dB • DL SINR above −7 dB when all HAPSs transmit signals at maximum power (full buffer). These SNR constraints are defined based on the required SNR, including a 15 dB margin. While ensuring the area coverage by the SNR constraints, the SINR constraint for DL is defined to avoid degrading communication quality because DL traffic is dominant over UL traffic. Here, we use coordinated multi-point (CoMP) reception for UL, which is commonly used in mobile wireless communications systems, e.g., 4G LTE. The UL CoMP allows us to obtain the diversity gain by receiving signals from UE in several neighboring cells, thereby improving the UL SNR, particularly at cell borders, which extends the coverage area. In this paper, we assume UL CoMP with two cells.
The SINR and SNR calculations are described in the following. The DL received signal level for the u-th UE from the c-th cell in dB is calculated as follows: where P t is a BS transmission power, and G u,c is the BS antenna gain of the c-th cell for the u-th UE, which is based on (2). The G UE is the antenna gain for UE, and L u is the free-space path loss between the u-th UE and HAPS, which is calculated as: where r u is the distance between the u-th UE and HAPS, and λ is the wavelength. Here, N values are calculated for the u-th UE, where N is the number of cells in a HAPS. The cell to which the u-th UE connects is defined as the cell with the highest received level among N cells. Let c u be a cell index to which the u-th UE connects. The total interference power of the u-th UE, which is connected to the c u -th cell in dB, is calculated by as follows: where β ∈ [0, 1] is an activation factor representing how many RBs are used in adjacent cells. We assume an activation factor of 50% (β = 0.5), which means that 50% of resource blocks are utilized in neighbor cells for the interference calculation. The DL SNR and DL SINR of the u-th UE, which is connected to the c u -th cell in dB, are γ SNR,u = 10 log 10

  10
Ru,c u 10 where σ 2 is the noise power. The noise power in dB is calculated as where B is the bandwidth, and F is the noise figure. Similarly, the UL received signal level of the u-th UE in the c-th cell is calculated as follows.
where P t is the UE transmission power. The cell to which a UE connects is determined by the DL received level, i.e., the cell with the maximum DL received power, which denoted c u . When the second strongest signal power from the u-th UE is denoted R u,c u where c u = c u , the UL SNR is calculated assuming UL CoMP as follows:

III. OVERVIEW OF DYNAMIC CELL OPTIMIZATION
Dynamic cell optimization attempts to optimize individual cell antenna parameters for a multi-cell configuration in terms of communication quality and capacity considering the user distribution. Here, we assume a cylindrical phased array antenna [5], wherein the horizontal and vertical beamwidths and directions can be controlled independently for each cell. The user distributions can be estimated by the signal strength obtained at the BS, the RTT, and the direction of arrival estimation using an UL signal. In addition, GPS feedback from users is another solution. Note that, we assume ideal estimation of the user distributions.

A. CELL OPTIMIZATION PROBLEM UNDER NON-UNIFORM USER DISTRIBUTION
The cell configuration shown in Fig. 2 assumes a uniform user distribution, and cells in the same layer have the same beamwidth and vertical direction angle, where the horizontal direction angles are spaced equally. In contrast, each cell should be controlled using different values for the same parameters compared to other cells while considering non-uniform user distributions. As described in Section II, each cell has four parameters (Fig. 3). Thus, the total number of parameters is four multiplied by the number of cells, e.g., the number of parameters is 36 for nine cells and 84 for 21 cells. Note that an increasing number of parameters results in a huge search space; thus, parameter optimization becomes difficult in some cases, even when using a GA or an exhaustive search technique. Thus, to address this problem, we propose a two-step dynamic cell optimization algorithm for non-uniform user distributions. The proposed method comprises both search space reduction and antenna parameters optimization steps. The proposed method reduces the search space using a co-evolutionary algorithm to reduce the number of parameters to optimize in a single GA execution and focuses on a promising search range for each parameter in the search space reduction step.

B. CO-EVOLUTIONARY DYNAMIC CELL OPTIMIZATION
If the distance between cells A and B is large, changing the cell A's parameters has a little impact on cell B. In other words, cell parameters with a large distance between the cells exhibit low dependency on each other. Thus, a method that divides multi-cells into several groups (i.e., sub-areas) and optimizes each sub-area in an autonomous and distributed manner is possible. However, an autonomous and distributed approach is insufficient because it cannot consider the dependence among sub-areas and coverage conditions as an entire area. To address these problems, we propose a cell optimization method with co-evolution, where each sub-area is optimized iteratively [21]. This method can reduce the number of parameters to optimize by dividing multi-cells and enables both local optimization in a sub-area and overall optimization by going around the sub-areas. This method reduces the number of combinations compared to a technique that does not divide multi-cells. Fig. 4(b) shows an example of multi-cell division, where nine cells are divided into three groups. When we focus on a given sub-area, only the corresponding parameters are optimized. In Fig. 4(b), the number of parameters is reduced from 36 to 12, as shown in Fig. 5. In another example, only four parameters are optimized for nine sub-areas. In this manner, 4N parameters are reduced to 4(N /M ) parameters according to the number of sub-areas M (note that N is divisible by M ). The proposed method optimizes the parameters based on the division of multi-cells as follows. This process continues until all M sub-areas in the area are selected. In addition, this process can repeat T times such that steps 1-3 are iterated MT times in total, where T is the number of optimization cycles.
Sub-area division is also effective in terms of reducing the calculation time of the objective function in the optimization process. As described in Section II-C, the received signal power is calculated for each cell per UE. The received signal power related to non-target sub-areas can hold the same value throughout the GA process because the parameters related to those sub-areas are constant. Although the actual effect of calculation time reduction for an objective function depends on its implementation, the required calculation time is reduced as the number of sub-areas increases.

C. SEARCH SPACE REDUCTION BASED ON MARGINAL HISTOGRAM
In addition to co-evolution-based sub-area division, focusing on a promising search range of each parameter is effective to reduce the combinations of parameters. However, it is difficult to predefine search ranges for various user distributions; thus, we also propose a method to reduce the search ranges for each parameter based on marginal histograms generated from the population optimized with a large step size. The basic concept of the proposed method is that the candidate solution distributes in a promising range because the GA evolves the population, and the characteristics of the population as a whole improve with each generation. This method is partly inspired by a probabilistic model-building genetic algorithm (PMBGA) [22] or estimation of distribution (EDA) [23] in which an offspring is generated based on a population probability distribution. For example, a previous study [24] utilized a marginal histogram in the context of PMBGAs. Note that the proposed method differs from PMBGAs in that it produces the probability distribution of each parameter with a large step size optimization to reduce the search ranges.
Prior to applying the search space reduction algorithm, each sub-area is optimized with a larger step size than the optimization step. In the search space reduction step, the optimization does not necessarily have to converge because we only need to know the probability distribution of each parameter. Thus, the number of generations for search space reduction can take an arbitrary value, which affects the accuracy of the probability distribution in the population.
The basic procedure of the proposed search space reduction algorithm is described as follows.
1) A marginal histogram about the given parameter is generated and find an index s of the peak value.  2) Extend the range by ±1 and calculate the sum of probability p sum in the range. 3) Repeat step 2 until the sum of the probability p sum exceeds a threshold t. 4) Define a new range based on the range obtained in step 3. The above procedures are illustrated in Fig. 6, which shows an example of reducing the search range where the original range is −60 to 60. In addition, a flowchart of the proposed search space reduction algorithm is given in Fig. 7. In step 1, a marginal histogram is generated by focusing on a single parameter. Here, each bin in the histogram corresponds to the candidate value for each parameter. Generally, considering a joint probability of p(X 1 , X 2 , . . . , X 4N ), the marginal probability of X 1 is represented as follows: The probability of the i-th bin in the histogram is represented as: where c i and N pop are the number of elements of the i-th bin and that of populations in a GA, respectively. The sum of the probability is calculated as where m is the number of iterations for step 2. Note that when s + i becomes less than one or exceeds the maximum index, we assume p s+i = 0 for the calculation of (14). If probability p s is greater than the given threshold in step 2, we skip steps 2 and 3, and a new range is defined based on a range from s−1 to s + 1. The redefined range retains the original lower/upper range if the peak index s is one or the maximum in step 2.
Although the threshold can take a value up to 1, the effect of reducing the search space reduction becomes small with a high threshold. Here, a small threshold reduces the search space excessively, and, as a result, the optimal solution will not be found. As shown in Fig. 6, when the threshold is set to 0.68, s = 8, and m = 2, the new range is from −16 to 16.  space reduction step (step 1), where the search range for each parameter is optimized with a large step size. The search space reduction algorithm described in Section III-C is applied only when T = 1. The rest of the optimization cycles (T ≥ 2) are dedicated to the optimization step (step 2), where the reduced search ranges for each parameter in step 1 are optimized with the desired step size. Each sub-area is co-evolutionary optimized regardless of the step as described in Section III-B. The algorithm starts with the initial parameters which are optimized for a uniform distribution. For example, the initial parameters for N = 9 with a radius of 50 km are shown in Table 2, which corresponds to Fig. 2 (a). The algorithm terminates when the optimization cycle T exceeds the predefined optimization cycles T max .

IV. GENETIC ALGORITHM AND PROBLEM FORMULATION A. OPTIMIZATION BY GA IN THE PROPOSED METHOD
In the previous subsection, the overall process of the proposed method was described. In the optimization step, we employ the GA to optimize parameters because GAs can be parallelized easily, thereby improving optimization speed. In the GA, a population of candidate solutions called individuals goes through operations inspired by evolutionary processes, e.g., crossover and mutation, to realize a better solution.   Fig. 4 (b)).
In the first generation, the parameters for a uniform user distribution are incorporated into the population to reduce optimization time. Here, the population size is set to 100. For each individual in the population, the fitness value is calculated by an objective function. We define the objective function as the sum of the square root of the throughput for all users as an example [25]: where N UE , B, N c u , and γ SINR,u represent the total number of users, the bandwidth, the number of users in the cell where the u-th user connects to, and the DL SINR of the u-th user, respectively. We select the objective function in (15) because, if the sum throughput is used as the objective function, some users with large throughput may affect the objective function significantly. After calculating the fitness values based on (15) for all individuals in the population, the best 5% of the population is preserved as elite offspring, 80% of the new population is produced via a crossover operation, and 15% of the new population is produced by a mutation operation. Here, implementation of the crossover and mutation operations is based on the literature [26]. The objective function is modified into a penalty function in a manner that considers a violation of constraint functions based on the literature [27].
The coverage conditions (Section II-C) are evaluated under the assumption of uniform distribution such that the coverage conditions are satisfied for areas without users.
The GA continues producing new generations to obtain a better solution, and termination occurs when the stopping criteria are satisfied, i.e., the best penalty value remains unchanged for certain generations, which we refer to this period as stall generations.

B. PROBLEM FORMULATION
Here, we describe the problem formulation. Individuals are defined as a combination of antenna parameters to optimize. A vector for an individual is represented as p = [p 1 , p 2 , · · · , p 4N ] where the size of the vector is determined by the number of parameters. The k-th element of p is a positive integer value, and its maximum number is the number of candidates with corresponding parameters. For example, when a horizontal direction that ranges from −60 to 60 degrees with a step size of 2 degrees is optimized, p k takes a value of 1 to 61, where the index of 1 corresponds to −60 degrees. The default range for each parameter is given in Table 3. The parameter n is 3 for the inner layer and 6 for the outer layer for nine cells, as shown in Fig. 4. In reference to the literature [6], we assume a step size of 2 degrees in this paper.
When the total number of combinations is represented as J , the pool of available parameters is denoted P : {p 1 , p 2 , · · · , p J }. The optimization problem is to find the best parameter combination P out of P considering the coverage conditions described in Section II-C, which is formally expressed as follows: where DLSNR 1% , ULSNR 1% and DLSINR 1% are the 1st percentile values of the DL SNR, UL SNR, and DL SINR of users in the area. Note that the DL SINR is calculated assuming a full buffer, and the size of p is reduced to 4(N /M ) when the proposed method is applied.

A. EVALUATION CONDITIONS
In this evaluation, we assume nine-cell and 21-cell configuration as shown in Fig. 2. Here, the area radius is 50 km with a single HAPS. The number of users is 3000, which are connected to the cell with the highest receiving power. The primary system parameters are listed in Table 1 (Section II).
In addition, important variables in this paper are listed in Table 4. The actual user distribution varies with the type of city and landform, e.g., mountains and rivers. Thus, it is crucial to evaluate the proposed method with a variety of user distributions. Here, we use population data provided by the AGOOP Corp. [28] to generate the user distribution for our VOLUME 10, 2022  simulation. The population data counts the population size in a mesh of approximately one square kilometer. Here, the center coordinates are selected randomly, and 3000 points are selected randomly from all points in the area with the selected center coordinates. Although there may be millions of users in the given area, not all users necessarily need to be considered to calculate the objective function in (15) because the total throughput is unaffected by the number of users due to the factor B/N c u in (15) unless the number of users is not too small. We set the number of plots to 3000 in the evaluation to reduce simulation time under the condition that the center of gravity of the extracted plots from the original distribution retains that of the original distribution. We consider 100 different user patterns in this study. Fig. 9 shows examples of non-uniform user distributions, where the map is based on the maps published by the Geospatial Information Authority of Japan.

B. RESULTS
As stated previously, we evaluated the proposed method for nine-cell (N = 9) and 21-cell (N = 21) configurations compared to a method that does not apply sub-area division (M = 1) as a benchmark. Here, the results are evaluated in terms of the average improvement of the objective function value, where the reference value is calculated with the parameters for a uniform user distribution (Fig. 2). The number of average performed combinations in the optimization is also evaluated. In the proposed method, the GA is executed M times larger than in the method without sub-area division; however a complexity evaluation according to the number of performed combinations is good because no additional computations are required for sub-area division. In contrast, the complexity of the search space reduction algorithm is O(n 2 ) because it depends on the number of parameters and the step size in the search range. However, in our simulation, the computational time of search space reduction algorithm is negligible compared to that of the calculation of the objective function. Thus, we consider the average number of performed combinations even in the evaluation of step 2.
First, it is necessary to assess the effect of optimizing cell configurations in terms of the objective function (15) according to non-uniform user distributions. Fig. 10 shows examples of cell configuration optimized for non-uniform user distributions, where a dark yellow areas have a strong received signal level. Table 5 shows the parameters corresponding to case shown in Fig. 10(a). The results show that some cells with intense received signal power are modified to cover user clusters. To characterize the generated user distributions, we define the ratio of the number of meshes with users to that of all meshes in the area as a population distribution ratio. Fig. 11 shows the relationship between the population distribution ratio and the improvement rate of the objective function, where the reference value (0%) is calculated with the parameters for uniform user distributions. As an example of N = 9, the results of M = 1 with T = 3 and 500 stall generations, M = 3 with T = 3 and 200 stall generations, and

TABLE 5.
Parameters for a non-uniform user distribution for N = 9 ( Fig.10 (a)). M = 9 with T = 3 and 75 stall generations, are plotted. These result demonstrate that a lower population distribution ratio results in a higher probability of a high average improvement rate. Users tend to be crowded in a specific area in the user distribution with a low population distribution ratio due to mountains and other landforms. Thus, we conclude that dynamic cell optimization is effective in non-uniform user distributions, particularly with a low population distribution ratio.
Next, we present the results of the proposed method without search space reduction (step 1) to observe the effect of sub-area division. Fig. 12 shows the average improvement rate of the objective function for N = 9 with M = 1 as a benchmark, where all parameters are optimized without sub-area division with different stall generations settings. As can be seen, when the stall generation is set to 50, the average improvement rate is approximately 5%, which is relatively low compared to that obtained with other stall generation settings. Thus, we consider that a stall generation  of 50 is insufficient for the search space. In some cases, the optimization fails to find a better solution than the initial parameters with small stall generation settings. By increasing the stall generation value, the search time becomes long such that the average improvement becomes approximately 17% with a stall generation value of 500. However, a large stall generation value, (>500) is not necessarily effective because the average improvement rate gradually converge to ∼17% by increasing the stall generation value in Fig. 12. Figs. 13 and 14 show the average improvement rate of the objective function for the proposed method where the number of sub-areas is set to M = 3 and 9 with different stall generations settings. As can be seen, the average improvement increases as T increases in both cases. In addition, the average improvement rate increases with increasing stall generations, where the maximum average improvement rates are approximately 21% for M = 3 with T = 5 and 16% for M = 9 with T = 5. The difference in the maximum average improvement rate between M = 3 and M = 9 is approximately 5%. In the M = 9, it is difficult to adjust the parameters with other cells because only a single cell is optimized at a time. Thus, the M = 9 case tends to fall into a local optimal solution. VOLUME 10, 2022     be found. For example, Fig. 13 shows that stall generation values of 100, 150, and 200 result in approximately the same average improvement rate for M = 3. A large stall generation setting may work as overhead because the GA continues to search for the period of stall generations when the GA finds a candidate solution. It is obvious that the average number of combinations for M = 9 in Fig. 17 is much less than that of M = 3 in Fig. 16. The number of parameters optimized at a time is 12 and four for M = 3 and 9, respectively. This makes the search space for M = 9 much smaller than that for M = 3 because the number of combinations grow exponentially according to the number of parameters. However, the M = 9 case appear to drop in local optima due to fast convergence.
To summarize the above results, Fig. 18 shows the relationship between the average improvement rate and average combinations for M = 1, 3, and 9 with different numbers of T . Here, each point corresponds to the stall generations settings in Figs. [12][13][14] where the leftmost points correspond to the minimum number of stall generations for each case. The more a plot goes to the lower-right area, the better the performance where a high average improvement  rate is achieved with a low average number of performed combination. The M = 3 with T = 1 case achieves an average improvement rate of 17% with approximately 65,000 combinations, and the M = 1 with T = 3 case achieves the same average improvement rate with approximately 375,000 combinations. This demonstrates the effectiveness of the proposed method in terms of reducing the search space. In addition, the maximum average improvement rate obtained with M = 3 is better than that of M = 1 by 3.36% with T = 3. This is because the better solution becomes easier to find by decreasing the search space when M = 3. Although M = 9 is also effective in terms of reducing the number of combinations, the average improvement rate is limited. The point located in the lower-right area most is the M = 3 with T = 2 case. Thus, M = 3 is effective in terms of both average improvement rate and average combinations compared to the M = 1 and M = 9 cases.
We also present the results of N = 21 for a different type of sub-area division. Fig. 19 shows two types of seven sub-areas for N = 21, and Fig. 20 shows the relationship between average improvement rate and average combinations for M = 1, 3, 7, and 21 with different T values. Each point in Fig. 20 corresponds to the stall generations settings in Table 6, where the leftmost points correspond to the minimum number of stall generations for each case. All types of sub-area divisions worked to increase the average improvement rate of the objective function compared to the case without sub-area division (M = 1). The M = 7 case is the best in terms of the average improvement rate of the objective function among different sub-area division settings, where the lines of two types of seven sub-area divisions nearly overlap. The average improvement rate increases by approximately 16% with 300,000 performed combinations by applying sub-area division with seven groups compared to the M = 1 case. For 21 cells without sub-area division, the number of parameters is 84, which makes this case even more difficult to optimize than the nine cell case. The M = 21 case reduces the average number of combinations; however, a decrease in average improvement of approximately 3.5% compared to the M = 7 case is observed. In addition, the M = 3 case exhibits a lower average improvement rate than the M = 7 and M = 21 cases. It is considered that the optimization is insufficient because the number of parameters remains large.
Next, we evaluate the proposed method with search space reduction (step 1) for the N = 9 with M = 3 and N = 21 with M = 7(1) cases, which are the best sub-area division settings in terms of average improvement rate of the objective function. As a benchmark, the results without step 1 are considered to assess the effect of step 1. Here, all parameters are optimized with a step size of 4 degrees in the search space reduction step. Fig. 21 shows the relationship between the average improvement rate and average number of combinations for the N = 9 with M = 3 case for different settings for the number of generations in the search space reduction step. Here, the threshold for the search space reduction algorithm is set to 0.68. Each plot with search space reduction in Fig. 21 corresponds to the stall generation settings of 25, 50, and 100. The number of generations for the search space reduction step is 25, 50, and 100, where optimization with the GA in the search space reduction step stops after this period. As can be seen, search space reduction with the number of generations for the search space reduction step of 100 shows nearly the same average improvement rate compared to the case without search space reduction while also reducing the average number of combinations. For example, the average number of combinations is reduced by 43% with search space reduction to achieve a 20.2% average improvement rate. Search space reduction with 25 and 50 generations for the search space reduction step exhibits a reduction in the average improvement. We assume that a small number of generations in the search space reduction step is insufficient to produce effective histograms to truncate unnecessary search ranges. Note that the average improvement rate for T = 1 is relatively low with search space reduction because the optimization is performed with a step size of 4 degrees and stops before convergence. Fig. 22 shows the relationship between an average improvement rate and an the average number of combinations for the N = 9 with M = 3 case with different threshold settings. Here, the number of generations for search space reduction is 100. The result demonstrates that the average improvement rate increases as the threshold for step 1 increases although the average number of combinations does not exhibit a clear difference among the different threshold values. We assume that setting a small threshold reduces the search space of each parameter excessively, thereby resulting in missing a good solution and reducing the average improvement rate.
Finally, we discuss the results obtained by the proposed method with step 1 for 21 cells. Fig. 23 shows the relationship    the average number of combinations approximately the same as the case without step 1 (Fig. 23). Similar to the results obtained with N = 9 in Fig. 22, the average improvement rate increases as the threshold for step 1 increases in Fig. 24. However, the proposed method with step 1 for a threshold of 0.9 does not demonstrate an obvious reduction to the average number of combinations compared to the case without step 1.
In summary, sub-area division is effective for the average improvement rate by reducing the number of parameters for both the nine-cell and 21-cell cases. However, we find that the search space reduction algorithm is only effective for nine cell configurations. A versatile search range reduction algorithm that can also be effectively applied to 21-cell configuration will be future work.

VI. CONCLUSION
In this paper, we have proposed a two-step dynamic cell optimization algorithm to maximize the objective function defined as the sum of the square root of throughput for non-uniform user distributions. User distributions differ depending on location; thus, it is beneficial to optimize the antenna parameters of cells accordingly. In addition, optimizing cell configurations for non-uniform user distributions cause exponential expansion of the search space, thereby increasing the difficulty of optimization even when using a GA in some cases. The proposed method comprises both search space reduction and optimization steps. The proposed method divides multi-cells into sub-areas comprising neighboring cells, which reduces the number of parameters that must be optimized. In the search space reduction step, the search range for each parameter is reduced based on a marginal histogram obtained via optimization with a large step size. Our simulation results demonstrate that cell configurations optimized for non-uniform user distributions outperform those for uniform user distributions in terms of the objective function value. The results also demonstrate that sub-area division is effective for nine-cell and 21-cell configurations, and that the search space reduction step can reduce the number of performed combinations for nine-cell configurations.