Foam Evolution Inspired Modeling for Staged Construction of Ultra-Dense Small Cell Networks

Small cells (SCs) are expected to be ultra-densely deployed in or close to the traffic hot-spots in the fifth generation (5G) mobile networks to provide wireless capacity cost-effectively. Traffic hot-spots change over time, which means SCs cannot be deployed in a one-off manner as macrocells normally do, rather they should be constructed in a staged process. Hence, mathematical models that capture the time-varying staged-construction process, are urgently needed for operators to effectively predict the construction period, but are currently lacking. In this paper, inspired by the foam bursting process-a natural phenomenon that can be observed in daily life such as hand-washing, we first propose a novel model that can predict the time-varying expectation and logarithmic variance of SC coverage areas. Then, we verify the model by real network deployment cases. Additionally, in order to extract parameters from historical base station deployment data, a parameter estimation algorithm is designed and verified. The findings of the paper reveal that mobile operators should construct ultra-dense SC networks in a staged manner like how larger foams split into smaller ones.

network to the public after completing the initial construction stage. Therefore, the operating revenue obtained from the earlier incomplete-network opening can cover parts of the network capital expenditures (CAPEX) and the operating expenditures (OPEX) [11] of the subsequent construction stages.
We define this new construction pattern as a staged construction (SCON) process, in which the expectation and variance of the BS coverage area change as the construction progresses with stages. This temporal change of BS coverage area distribution is closely related to the network performance. Therefore, mathematically modeling the variation of the time-varying BS coverage area distribution is crucial for operators to effectively predict the construction period required for network construction using the SCON approach. Specifically, this mathematical model can help operators formulate the corresponding SCON plan to finally achieve the target network performance.
Existing studies on dynamic networks have been proposed in the field of self-organizing networks and BS on/off switch handling [13]- [15]. In these researches, to improve the ultra-dense small cell network's energy efficiency, researchers have proposed different mathematical models that can automatically switch the state of SCs to respond to the fluctuating network traffic. These mathematical models can improve the network performance, such as network energy efficiency or the cell-edge throughput for the fully constructed networks. However, existing models ignore the dynamic network changes in the construction process [13]. To the best of our knowledge, the mathematical model to describe the under-constructed network, i.e., a network with time-varying BS distribution, has never been investigated before.
Interdisciplinary research often plays an essential role in accelerating disciplinary discovery in blank study fields such as physics-based or biology-based engineering [16]. As the first SCON network modeling, in the absence of existing research on the BS construction process mechanism, we introduce the process with similar evolution from physics-based engineering, i.e., the inverse foam bursting process, as a reference for the BS construction process [17]. The inverse process of bubble change during the foam bursting process is highly similar to the SC coverage change in SCON. During the foam bursting process, the bubble evolves from initial small average size and uniform distribution to eventual large average size and non-uniform distribution. By contrast, during the SCON process, the coverage area of SCs evolves from initial large average value and non-uniform distribution to eventual small average value and uniform distribution.
Therefore, to fill the gap in the dynamic modeling of ultra-dense small cell deployment using the SCON approach, we propose a novel time-varying mathematical model of BS distribution inspired by the foam distribution model to lay a foundation for this research field. Moreover, we provide a case study for operators to effectively make a trade-off between the construction period and the service probability with the proposed model while designing a new SCON network.
The specific contributions of this work are summarized as follows: • The expectation and the logarithmic variance of BS coverage areas, which are modeled as variables of the time-varying coverage area distribution in a point process, are employed as metrics of temporal network performance.
• The inverse foam bursting process is applied to the BS SCON model. The spatial-temporal mathematical expressions of the average and the logarithmic variance of BS coverage areas are derived.
• A Kolmogorov-Smirnov statistic [18] based estimation algorithm is proposed to extract parameters of the proposed model from actual Long-Term Evolution (LTE) BS spatial-temporal distributions.
• A method for computing the service probability of a new under-construction network using the proposed model is presented. The rest of this paper is summarized as follows. Section II presents the framework of the time-varying ultra-dense SC network. Section III presents the modeling of the SCON network. Section IV presents the algorithms that can extract parameters from historical network construction data. Section V illustrates the validation based on the historical data and a case study for how operators weigh the construction period against the service probability with the proposed SCON model. Section VI concludes this paper.

II. TIME-VARYING ULTRA-DENSE SMALL-CELL NETWORKS
Existing mathematical models usually assume that the BSs' locations follow a Poisson point process (PPP) [19], [20]. Calculations of network performance given from these models are tractable. The area distribution of Voronoi cells in PPP was simulated mathematically in [21]. Moreover, a network performance framework was provided for this mathematical modeling [22]. Let V be the area of a generic PPP-Voronoi cell, then: where λ is the density of BSs, and c is the fitting constant. When c is 3.5, this distribution fits the 2D PPP. According to (1), the expectation (V textave ) and logarithmic variance (σ D[ln V ]) of V in PPP can be derived as: Equations (2) show a bijection between the expectation and variance of the BS coverage areas when BSs' locations follow PPP. However, positions of SC BSs are affected by complex factors such as the change of hot-spots and the constraints of the construction environments. Thus, the new hot spots are not uniformly distributed, which destroys the correlation between the expectation and variance of the BS coverage area. Therefore, closer to the actual BS distribution, we assume that the location of the BS follows a point process, where the distribution of Voronoi cell area V still follows (1), but c becomes a variable controlling the variance of V . It provides additional adjustability to the BS modeling without adding extra complexity to network performance estimations based on the parameters λ and c. In this paper, V [km 2 ] represents the BS coverage area of its Thiessen polygons [23] when Voronoi Tessellations is used, as shown in Fig. 1.
Based on an open-access database containing coordinates and completion time of LTE BSs in an urban region of China [24], we consider the Thiessen polygon area of each BS generated by Voronoi Tessellation as its ideal coverage area. We download an LTE BS database in SQLite file. The data includes LTE BS data of construction date and latitude-longitude coordinates and can be directly extracted by MATLAB. The construction period extracted from the actual data is 32 days, divided into 10 construction stages according to the practical constructed date of these BSs. We calculate statistics of the actual LTE BS coverage areas and their distribution at different construction stages. Fig. 2 maps BS coverage at both middle and late construction stages. At the unfinished stage of the construction process, operators give priority to cover the populated regions as they are not able to cover all areas in an urgent schedule. As a small number of BSs are expected to cover a large region and to serve massive users, both the expectation and logarithmic variance of the coverage area are large. By contrast, when the network construction is finished, to improve the service probability in all regions, more BSs are deployed. In that case, both the expectation and logarithmic variance of the BS coverage area become smaller when the distribution tends to be uniform.
Correspondingly, Fig. 3 shows that the PDF of coverage areas of real LTE BSs matches that of Voronoi cell area distribution in the modeled point process well. Note that the data of LTE BSs over a city as the empirical data is used as an empirical SCON data in this paper. These LTE BSs will be more clustered by the streets than the ultra-dense SC networks distributed in a smaller area, which will worsen the fitting between the actual data and the theoretical distribution. We foresee that the PDF of actual data will fit better with the theoretical PDF if the SCON data of the ultra-dense SC network can be obtained.
Also, Fig. 3 shows the following construction pattern: with the construction processing, the number of deployed LTE BSs increases, while both V ave (t) and σ (t) decrease. To quantitatively describe this pattern, we define V ave (t) and σ (t) as functions of time t. By deducing the formula of ∂V ave (t) ∂t and ∂σ (t) ∂t , we can derive the relationship between the distribution of BSs and time. A mathematical model that fits the trend of V ave (t) and σ (t) in actual SCON network data, can be applied to predict the BS layout at different construction stages in the new SCON model with the similar construction conditions, and thus to predict the network performance at these stages. Therefore, mathematically modeling the variation of the time-varying BS coverage area distribution is crucial for operators to effectively predict the construction period required, and thus operators can use this model to construct an SCON ultra-dense SC network.

III. FOAM BURSTING PROCESS INSPIRED SCON NETWORK MODEL
The mathematical model of the foam bursting process shows heuristic significance for us to simulate the time-varying distribution of BS coverage areas in an SCON process. Kostoglou et al. [17] studied the discipline of the time-varying expectation and logarithmic variance of bubble volume in the bursting process for a large amount of foam in a fixed volume. This foam bursting process is of foam evolving from high density and uniform distribution to low density and non-uniform distribution. The inverse process is enlightening to the coverage area distribution in SCON of ultra-dense SC networks. Besides, similar to the process of foam bursting, the variations of parameters in the BS construction process are discontinuous.
Note that the process of the foam bursting is not precisely the same as the process of SCON, as latter contains more complicated factors. The foam bursting process is caused by the physical mechanisms of the foams, while the BS construction process is caused by the construction planning and the change of hot spots within the construction scope. The foam bursting process is influenced by two main mechanisms. Firstly, the gas in the bubble diffuses through the liquid membrane from small bubbles (high pressure) to large bubbles (low pressure). The second is that the liquid membrane between adjacent bubbles becomes unstable and eventually breaks, causing the bubbles to merge [25]. These two physical mechanisms advance steadily. In a foam bursting process, these variables, although complex, can be computed by certain formulas [17]. However, the BS construction plan designed according to the reverse process of the foam bursting process is still not the same as the final project result since factors such as weather and holidays affect the construction interval, thus affecting the staged construction of BSs.
Nevertheless, given the gaps in the SCON modeling study, we introduce the inverse foam bursting process as a reference for the BS construction process and ignore the unpredictable factors to make the model more tractable. Generally, investigation on the inverse foam bursting process will still shed light on how SC BSs are deployed due to the similarity of both processes. Therefore, we adopt the inverse foam bursting process model to lay a foundation in this research field.
Based on the foam bursting model [17], we assume the dynamic time-varying functions of the expectation and logarithmic variance for the coverage area of the ultra-dense SC BSs as and It is noted that, when studying the foam evolution, parameters k and q are used to comprehensively describe the changing stability of liquid film between bubbles in the bursting process and thus vary with time. However, parameters k and q varying over time will bring more uncertainties in the construction plan, thus lead to a higher computational complexity without increasing the predicting accuracy. In this paper, we decompose k and q into two influencing factors for the BS construction process: k affects construction speed, and q affects construction uniformity. We assume that k and q are static to simplify the analysis, ensuring no additional complexity is involved. With values of k and q, which can be extracted from actual BS data, we propose a mathematical model that can fit the distribution change of actual BS coverage area data in Section III, so that the distribution change of the staged constructed BSs in a similar construction scenario can be predicted.

IV. PARAMETER EXTRACTION ALGORITHM
To simplify (3) and (4), we approximate the differential equations by difference equations as: where t is the sampling interval of the iteration. This mathematical model requires initial values for engineering prediction. The initial values we use are the expectation and logarithmic variance of coverage area obtained at the first sampling time in the actual construction data, V ave (1) and σ (1). According to (2), (5), and (6), temporal variables λ(t) and c(t) can be computed. In addition, λ(1) and c(1) derived by V ave (1) and σ (1) are also the initial values to simulate the expected theoretical cumulative distribution function (CDF) in the parameter estimation as follows.
According to Kolmogorov theorem [26], the larger the sample size is, the closer the distribution of actual data becomes as compared with the predicted distribution of the fitting model. Thus, on the premise that the actual data samples are reliable enough, the larger the sample size of the initial data is, the more accuracy of the parameter estimation will be. However, the initial values of V ave and σ only depend on the initial design of the network by operators in practical engineering applications. In this paper, we mainly study the change of BS deployments in the SCON process. For the initial low-density network deployment designing, many researchers have made relevant works, such as [13], [27]. Operators can refer to these works according to their needs.
To extract values of k and q from actual SCON network data, we design an algorithm to estimate these parameter values of k and q.
First of all, we design an algorithm to judge whether a set of parametersk andq ensure that the prediction model fits the actual BS construction data at each construction stage. The parametersk andq can be substituted into the difference equations (5) and (6), and we can obtain the theoretical V ave and σ at each construction stage, so as to the simulated distribution at each construction stage. The infinite norm between the theoretical CDF of simulated BS distribution, F V , and the CDF of actual BS distribution at each construction stage, F V , is calculated as the Kolmogorov-Smirnov statistic [18] to indicate the accuracy of parameters. The Kolmogorov-Smirnov statistic is denoted by F V −F V ∞ , which fits both V ave and σ in the estimation. According to Kolmogorov theorem [26], the square root of the sample size is inversely proportional to the distance between the empirical CDF and the CDF of the estimated model, so the square root of the sample size is used as the weight of the Kolmogorov-Smirnov statistic of each construction stage for estimation, denoted as D. The model with the smallest Kolmogorov-Smirnov statistic between simulated distribution and actual distribution in all construction stages is considered to have the best parametersk andq, which means the average of D at all construction stages, denoted as D mean ,

Algorithm 1 D mean Calculating Function
Require: • Sampling time N t ; • A set including all BS coverage areas d V and their CDF F V (N t ) at each sampling time; • V ave (N t ) and σ (N t ) at each sampling time; • The values of parameters k and q, denoted ask and q. Ensure: • The value of D mean that is derived fromk andq. 1: k ←k; 2: q ←q; 3:V ave ← V ave (1); 4:σ ← σ (1); 5: for n = 1 : length(N t ) do 6: d V at N t (n); 7: if n > 1 then 8: t ← N t (n) − N t (n − 1); 9: for m = 1 : N t (n) do 10  should be the smallest. The detailed steps of the function that computing D mean from parametersk andq are presented in Algorithm 1. When comparing the actual BS construction data with the fitted model, the CDFs of the fitted model and the corresponding actual CDFs at 10 construction stages are compared respectively.
The best estimation results of k and q from the potential range ofk andq are k and q that can minimize the value of D mean . To constrain the minimization of the D mean , denoted as D best , the interior-point approach [28] is applied to quickly compute the optimal solution in the estimation interval of parametersk andq, as shown in Algorithm 2.

A. HISTORICAL DATA-BASED VALIDATION
To verify the parameter estimation algorithm, we list the D mean results with different values of parameters k and q, as shown in Fig. 4. The minimum D mean = 1.85 is
• The initial values k 0 , and q 0 , ensuring that k 0 ∈k and q 0 ∈q. Ensure: • The minimum D mean that can be derived in the potential range, denoted as D best ; • The k and q that derive D best .  5: k ← x(1); 6: q ← x(2); 7: return k, q, D best found when k is 0.05 and q is 0.91, illustrating that the Algorithm 2 can effectively find the optimal solution to minimize D mean .
By substituting the values of the bestk andq into (5), (6), we have: where the path loss exponent is α = 4, and the target SINR isγ = 0 dB. Calculation of p service is referred to Equation (5) in [22].
Equations (7) and (8) can be used to estimate the time-varying expectation and logarithmic variance of BS coverage areas in similar construction scenarios. Then, we predict the trends of V ave and σ from the initial V ave as V ave (1) and initial σ as σ (1) by iterating over (7) and (8) with t as 1 day.

B. CASE STUDY
Combining (2), (7) and (8), the time-varying distribution of BS coverage areas can be derived. The actual BS construction data is thus converted into an empirical model. By combining this empirical SCON model with the existing solution for calculating the network service probability based on the Voronoi cell area distribution [22], we obtain the service probability of an under-construction interference-limited network for different MU density λ u , where the path loss exponent α = 4, and the target signal to interference and noise ratio (SINR) γ = 0 dB, as shown in Fig. 5. Note that the actual BS data is used to simulate the ultra-dense SC network, so that the densities of mobile users (MUs) are set to be 1 dB to 3 dB lower than the BS density (the BS density λ = 6.97 BSs/km 2 at the end of the construction when t = 100 days). When the operator set the target service probability (95% in Fig. 5), the expected period of a network construction scheme can be deduced through this model to evaluate whether the construction scheme, determined by k and q, is acceptable for this new SCON project. For example, if the target MU density of the network is 0.5 UEs/km 2 and the scheduled construction period is 80 days, the operators reject this construction scheme. On the contrary, if the target MU density of the network is 0.05 UEs/km 2 , and the scheduled construction period is 10 days, then this construction scheme can be considered acceptable.
It should be acknowledged that the actual BS density in the ultra-dense SC network will be much higher than the density of the actual BS data used in this paper. The total cost of every construction stage estimation for the SCON network and a one-off deployed network with the same construction scale. The CAPEX per cell is set as capex = $5000, and the OPEX per cell is set as opex = $1000. The target coverage range of the network S = 1km 2 . The actual cost data is referred from [11].
The proposed dynamic model for SCON SC network performance evaluation will evolve once the operators own the SCON data of an ultra-dense SC network. Fed with big SCON data, the proposed model lays a foundation to provide the design basis for the future ultra-dense SC network SCON projects.

C. COST ESTIMATION
In a network, according to [11], the CAPEX of each cell can be decomposed into the cost of the cell and the cost of the core network working on the cell, while the periodic OPEX per cell includes marketing costs, electricity bill, venue rental, back-haul rental, as well as the hardware and software maintenance costs. The cost estimation model can be simplified when assuming the CAPEX and OPEX of each cell as fixed values, denoted as capex and opex. Clearly that the CAPEX and OPEX of the under-constructed network at the construction stage t are related to the constructed BS number N (t) and the under-constructed BS number N (t + t) − N (t) at this construction stage: where N (t) = S × λ(t), and S refers to the target coverage range of the network. Combining (9), (10) with the actual cost data referred from [11], we estimate the total cost of every construction stage as the sum of CAPEX(t) and OPEX(t) of the empirical SCON model generated in Section V.B. The construction is assumed to be finished in day 100. The estimation result of cost per stage is shown in Fig. 6, with the costs of a one-off deployed network with the same construction scale plotted as comparisons. Fig. 6 clearly shows that a one-off deployment strategy results in significant initial investment requirements. In contrast, the SCON strategy would dilute this massive investment by a more extended construction period, reducing the operators financing difficulty. Moreover, all SCs have been deployed at the beginning of the one-off construction. That makes the OPEX of the one-off deployed network higher than that of the SCON network. The OPEX of the completed SCON network rises to the same as the one-off deployed network.

VI. CONCLUSION AND FUTURE WORK
In this paper, a mathematical model to describe the time-varying distribution of SC BS coverage areas in an SCON network has been proposed. The model was inspired by the foam bursting process. Moreover, parameters of this time-varying model have been estimated using actual BS data, and an empirical model has been obtained for future SC network deployment using the SCON approach. Combined with the existing network evaluation solutions, this mathematical model can help operators effectively evaluate whether a construction scheme generated by a completed SCON data is acceptable for designing new SCON networks. Since no available SCON data of ultra-dense SC networks exists at this moment, the proposed model is validated by comparing the distribution of BS coverage areas obtained by this mathematical model with actual LTE networks, which were built at a limited number of stages. With the large-scale deployment of ultra-dense SC networks, the increasing amount of reference data will make this model more accurate.
In the future, our model will evolve as the SCON data access increase. Thus, we will continue to study the SCON process mechanism in combination with economics, psychology, geography, and other interdisciplinary subjects. Moreover, since our work and the reviewed dynamic networks modeling in the field of self-organizing networks and BS on/off switch handling are focusing on the network at different construction states, we will combine our work with existing dynamic network models to optimize the network performance at different construction stages.