Joint Optimization of Interference Coordination Parameters and Base-Station Density for Energy-Efficient Heterogeneous Networks

Heterogeneous networks (HetNets), consisting of macro-cells and overlaying pico-cells, have been recognized as a promising paradigm to support the exponential growth of data traffic demands and high network energy efficiency (EE). However, for two-tier heterogeneous architecture deployment of HetNets, the inter-tier interference will be challenging. Time domain further-enhanced inter-cell interference coordination (FeICIC) proposed in 3GPP Release-11 becomes necessary to mitigate the inter-tier interference by applying low power almost blank subframe (ABS) scheme. Therefore, for HetNets deployment in reality, the pico-cell range expansion (CRE) bias, the power of ABS and the density of pico base stations (PBSs) are three important factors for the network EE improvement. Aiming to improve the network EE, the above three factors are jointly considered in this paper. In particular, we first derive the closed-form expression of the network EE as a function of pico CRE bias, power reduction factor of low power ABS and PBS density based on stochastic geometry model. Then, the approximate relationship between pico CRE bias and power reduction factor is deduced, followed by a linear search algorithm to get the near-optimal pico CRE bias and power reduction factor together at a given PBS density. Next, a linear search algorithm is further proposed to optimize PBS density based on fixed pico CRE bias and power reduction factor. Due to the fact that the above pico CRE bias and power reduction factor optimization and PBS density optimization are optimized separately, a heuristic algorithm is further proposed to optimize pico CRE bias, power reduction factor and PBS density jointly to achieve global network EE maximization. Numerical simulation results show that our proposed heuristic algorithm can significantly enhance the network EE while incurring low computational complexity.


Introduction
The exponential growth of data traffic demand, huge energy consumption and large amounts of global carbon dioxide emissions severely restrict the sustainable development of wireless cellular networks. According to the statistics, the data traffic volume demand in the fifth-generation (5G) wireless communication network will increase of 1000× by 2020. Moreover, the limited spectrum resources also constrain the network capacity improvement [1]. Therefore, the network energy efficiency (EE) which considers both spectral efficiency (SE) and energy consumption has been valued not only as an important network performance indicator for modern wireless networks, but also for the operational expenditure reduction and sustainable development [2].

Motivation
In HetNets consisting of MBSs and PBSs, as shown in Figure 1, the great difference of transmitter power leads to load imbalance between macro tier and pico tier. To address this issue, PBSs adopt cell range expansion (CRE) technology to enlarge the PBS coverage area by adding a positive bias on the reference signal received power (RSRP) of PBSs without increasing transmission power [5,6]. Unfortunately, CRE user equipments (UEs) located in the pico CRE region will suffer serious downlink interference from dominating MBSs, even causing the outage of control signals. As a result, it is important to mitigate the downlink interference to improve the wireless link quality between transmitter and receiver. Then the network EE can be improved. To mitigate the downlink inter-tier interference for HetNets, further enhanced inter-cell interference coordination (FeICIC) scheme has been specified in 3GPP Release 11 [7]. In FeICIC scheme, MBSs can transmit data to macro center UEs with low transmission power over certain subframes, termed as low power almost blank subframes (ABS), over which PBSs can schedule CRE UEs with reduced interference [8]. On the basis of FeICIC technique, for user association, the pico CRE bias will directly affect the value of user received RSRP from PBSs, which will eventually affect the number of CRE UEs associated to PBSs. As CRE UEs are scheduled in the subframes corresponding to the ABS of MBSs, the transmission power of MBSs in ABS, i.e., ABS power, will decide the interference degree suffered by CRE UEs. Therefore, pico CRE bias and ABS power are two important factors for improvement from different perspectives including resource management [19][20][21][22][23], FeICIC parameters optimization [24][25][26][27][28][29][30][31] and BS deployment strategy [32][33][34][35][36][37][38][39]. As for resource management, centralized resource allocation algorithms based on convex-optimization [19,20], graph-theory [21] or even game-theory [22,23] were proposed to achieve the maximum EE gain.

BS PBS
As for FeICIC parameter optimization, the ABS ratio dynamical optimization algorithm based on network load was proposed to enhance the network EE in [24]. For improving the network capacity and EE, pico CRE bias optimization problem was further developed combined with power control in [29]. Using the stochastic geometric approach, the expressions for SE and cell-edge throughputs have been derived as a function of the power reduction factor of low power ABS in [30]. To move one step further, pico CRE bias, ABS ratio and ABS power are jointly optimized by a robust EE optimization framework in [25]. In [31], it was proved that FeICIC can achieve a better gain in view of network EE, cell-edge throughput and user fairness compared with eICIC. In addition, the distributed algorithms based on the exact potential game framework for both eICIC and FeICIC optimizations were proposed to offer better network performance. The authors of [26] deducted and analyzed the EE coverage performance of FeICIC with extra pico CRE bias. In [28], FeICIC technique was applied to mitigate the inter-tier interference in HetNets deployed with unmanned aerial vehicles (UAVs). In such a scenario, the locations of UAVs, pico CRE bias and inter-cell interference coordination parameters were jointly optimized by using a genetic algorithm. In [27], the EE of HetNets with joint FeICIC and adaptive spectrum allocation was analyzed by the stochastic geometric approach. The research on FeICIC parameters optimization on the basis of stochastic geometry is relatively few in the existing literatures.
As for BS deployment strategy, general linear power consumption models were developed by means of linear regression in [38]. Meanwhile, the effects of MBS transmition power and on/off switching on instantaneous MBS power consumption were analyzed. A threshold of SBS density in ultra-dense HetNets was investigated in consideration of the backhaul network capacity and network EE in [32]. In [33], the authors came up with an approximation algorithm to solve the intractable user association problem by controlling the PBS density dynamically. The relationship between PBS density and network EE was analyzed under different UE density in [34]. It was proved that both PBS density and MBS density have a notable impact on the network EE in [35]. In [36,37], PBS density and MBS density were jointly optimized based on traffic-aware sleeping strategies and stochastic geometry to enhance the network EE, respectively. In [39], taking the traffic pattern variations into account, the BSs can not only adaptively switch on/off states but also can dynamically scale its transmit power according to network capacity demands. In this way, network energy consumption is reduced.
Despite the aforementioned research works, few works in the literatures focused on FeICIC parameters and PBS density joint optimization for network EE improvement in HetNets, which will be investigated and developed in this paper.

System Model
The traditional network models with a hexagonal grid cannot accurately match the actual network deployment. Under such deterministic grid models, Monte Carlo simulations consume huge amounts of time and resources to obtain the statistical results. Recently, a stochastic geometry model was proven to be a tractable analytical model for homogeneous networks and HetNets, where the location distribution of BSs is modeled as a spatial Poisson point process (PPP) [40]. Using PPP, the network performance, like signal to interference plus noise ratio (SINR) coverage [41], rate coverage [8], average rate [42], can be analyzed conveniently by theoretical derivation. Thus, we adopt a stochastic geometry model to model a two-tier HetNets consisting of MBSs and PBSs in this paper.
Let k ∈ {m, p} denote the subscripts of a tier, where k = m represents macro tier and k = p denotes pico tier. MBSs and PBSs are modeled as two identically independent distributions (iid.) PPPs Φ m and Φ p with density λ m and λ p in the Euclidean plane, respectively. The UEs are also distributed according to a different iid PPP Φ u with density λ u . The total spectrum bandwidth is defined as W.
To mitigate the downlink interference over the control channel from MBSs to CRE UEs, MBSs adopt the FeICIC scheme, where all the subframes are divided into protected subframes (PSF), i.e., low power ABS, and unprotected subframes (USF), i.e., non-ABS. The MBSs transmit data at reduced power ρP m on PSF, where 0 ≤ ρ < 1 is the power reduction factor. In fact MBS transmit power in USF and PBS transmit power will also have effects on network EE. However, to focus on the effects of PBS density, pico CRE bias and power reduction factor on network EE and also for analysis simplification, we assume that MBSs transmit power in USF and PBSs transmit power are set to be the maximum fixed power P m and P p , respectively. Let θ to be PSF ratio, which is defined as the proportion between the number of PSF subframes and that of all subframes. Each user is associated with the strongest BS according to the biased received reference signal received power (RSRP) at the user. In this paper, the association bias for MBS is assumed to be unity (B m = 1 = 0 dB) and that of PBS is pico CRE bias depicted as B p , where B p ≥ 0 dB.
Based on the user association strategy, all UEs can be divided into four different types: the type of PSF macro-cell UEs (MUEs) contains the users within the macro-cell center region, the type of USF MUEs includes the users outside the macro-cell center region, the type of PSF pico-cell UEs (PUEs) correspond the users located in the pico CRE region and the type of USF PUEs comprises the users scattered in the original coverage of pico-cell. As shown in Figure 1, we adopt the index l ∈ L = {pm, um, pp, up} to denote the indication of the above four types of users, respectively, where pm represents PSF MUEs, um denotes USF MUEs, pp signifies PSF PUEs, and up indicates USF PUEs. In HetNets with FeICIC, the UEs scheduling strategy for MBS and PBS is shown in Figure 2 According to the Slivyak theorem, there is no difference in properties observed either at a point of the PPP or at an arbitrary point [8]. Therefore, we can simply analyze a typical UE located at the origin. The received signal power of a typical user l from a BS of k tier at a distance of r k can be represented as P k hr −α k , where P k is the full transmission power of BS in k tier, h represents the channel fast fading gain, which is defined as Rayleigh distributed with average unit power, i.e., h ∼ exp(1). The term α denotes the large-scale path loss exponent, which is assumed to be the same in both macro tier and pico tier for convenient analysis. Hence, the SINR of a typical user l based on its user type can be depicted as: where I m and I p denote the full power aggregate interference from the macro tier and pico tier to UE l, respectively, σ 2 represents the thermal noise, ρ is the power reduction factor of MBS transmit power in PSF. P m and P p are the full transmission power of MBSs and PBSs, respectively. r m and r p are the nearest distance from MBSs and PBSs to a typical UE l, respectively. The notations used in this paper are shown in Table 1. Probability of a typical UE belongs to the user type l k c , k e , k p Factor of macro-cell center region, pico CRE region and pico-cell original coverage region r l Distance from a typical UE l and its serving BS f l (r) PDF of the distance between a UE and its serving BS λ l Density of BSs associated with user type l R l Mean achievable downlink data rate of a typical UE l W l Spectrum bandwidth allocated to a typical UE l N l , N total l Mean number of UEs with user type l in a Voronoi cell, total number of UEs with user type l P m,s , P p,s Static power of MBS and PBŜ P m Proportion between maximum transmission power of MBS and that of PBŜ P p Proportion between maximum transmission power of PBS and that of MBS λ p,m Proportion between PBS density and MBS density R total Total network throughput P total Total network power consumption ρ Approximate value of power reduction factor B * p ,ρ * ,λ * p Near-optimal value of pico CRE bias, power reduction factor and PBS density B * p , ρ * , λ * p , EE * Optimal value of pico CRE bias, power reduction factor, PBS density and EE

User Type Probability
Normally, the user type of a typical UE l can be decided by the relationship between the biased RSRP from its nearest MBS and its nearest PBS as follow: where B p is the pico CRE bias. In Equation (2), the conditions for determining user type can be further translated from biased RSRP based inequation into distance based inequation, which is shown as Equation (3).
where k c = B p P p / (P m ρ) 1/α , k e = B p P p /P m 1/α and k p = P p /P m 1/α are defined as macro-cell center region factor, pico CRE region factor and pico-cell original coverage region factor, respectively [8], determining the coverage bound of the macro-cell center region, the pico CRE coverage region and the pico-cell original coverage region, respectively. In order to obtain the probabilities of four user types, the following lemma is proposed.
Due to the locations of MBSs and PBSs follow two iid. PPPs, given two arbitrary coefficient values of n a and n b , the probability of n a r m < r p < n b r m can be expressed as: where λ m and λ p are the MBS density and PBS density.

Proof. The proof of Lemma 1 is presented in Appendix A.
Then, the probabilities of the PSF MUEs and the USF PUEs can be calculated similarly based on Lemma 1 and can be expressed by Prob r p > n a r m = λ m λ m +n 2 a λ p and Prob , respectively. Combining Equation (4) with the user association strategy in Equation (3), the probability of this typical UE belonging to the user type l can be defined as A l = Prob(l ∈ L), which is expressed as: In particular, if ρ = 0, then k c = ∞ and A pm = 0, which means that the PSFs are zero power ABS. The association probabilities of USF MUEs and PSF MUEs versus ρ with different B p are simulated according to Equation (5) in Figure 3. As shown in Figure 3, with the increase of power reduction factor ρ, the transmission power of MBSs over PSF will increase, resulting in the association probability of USF MUEs decreasing and that of PSF MUEs increasing. It can also be found that the sum of A pm and A um will decrease with the growth of B p , because more UEs will be offloaded into pico-cells with larger B p .

Distribution of Serving BS Distance
After a typical UE type is classified according to the user association strategy, the probability density function (PDF) of distance r between this typical UE and its serving BS can be obtained according to Lemma 2 as below.

Lemma 2.
On the basis of user association probability deduced in Equation (5), the probability density function (PDF) f l (r) of the distance r between the typical UE l and its serving BS can be derived as: Proof. The proof is shown in Appendix B.

Derivation of Energy Efficiency Expression
This section introduces our main analysis model and derives the closed-form expressions of the average achievable downlink rate, the network power consumption and the network EE, respectively.

The Ratio of PSF
PSF ratio can be denoted to be the proportion between the association probability of PSF PUE and the sum of the association probability of PSF PUE and that of USF PUE, as shown in Equation (7).
where the expressions of A pp and A up are obtained according to Equation (5).
Referring to Equation (7), the PSF ratio versus B p with different PBS densities λ p is depicted in Figure 4. As shown in Figure 4, with B p increasing, the pico CRE area will be enlarged. As a result, the PSF ratio will rise. Moreover, with the PBS density λ p increasing, the distance between PBSs will be smaller, which will limit the further expansion of pico CRE area, so that the effect of B p on θ will be weakened with λ p increasing.

Average Achievable Downlink Rate
Assume that the network system adopts full buffer model and the frequency resource is allocated to all UEs in the coverage of a BS equally. Thus, the mean achievable downlink data rate of a typical UE l can be denoted as: where W l is the spectrum bandwidth allocated to UE l. Specifically, when l ∈ {pm, pp}, W l = θW and when l ∈ {um, up}, W l = (1 − θ) W. N l is the mean number of serving UEs with user type l in a Voronoi cell and its expectation is According to the analysis above, we get Lemma 3 as follow: Lemma 3. The average achievable downlink rate of a typical UE l can be further represented by Proof. The proof is presented in Appendix C.

Corollary 1.
To further simplify the analysis, we ignore the noise, i.e., σ 2 = 0 and set the large-scale path loss exponent α= 4. In that case, corresponding to the user type, the average achievable downlink rate of four user types can be expressed as: where Q (τ, Proof. Set α = 4 and σ 2 = 0, then we can get ϕ l = 0 and Z (τ, Combining with Equation (6), the desired results in Equation (9) can be obtained.

Network Power Consumption
Generally, the BS power consumption comprises static power consumption and transmit power consumption [35]. The static power consumption is caused by signal processing, battery backup, as well as site cooling, and independent with the BS transmit power consumption. The transmit power consumption is determined by the transmission power of this BS and the load-dependent power consumption coefficient of this BS which is denoted as the number of its serving UEs. Define P m,s and P p,s are the static power consumption of each MBS and each PBS, respectively. With FeICIC scheme, the power consumptions of each MBS in PSF and USF can be expressed as P PFS m = P m,s + N pm ρP m and P UFS m = P m,s + N um P m , respectively. Similarly, the power consumptions of each PBS in PSF and USF can be given as P PFS p = P p,s + N pp P p and P UFS p = P p,s + N up P p , respectively. In PFS, the unit area mean power consumption P PFS can be expressed as: Similarly, the unit area mean power consumption in UFS P UFS can be obtained as follow: Hence, the network power consumption can be expressed as:

Network Energy Efficiency
The network EE can be defined as the ratio of the achievable network throughput to the network power consumption [37]. For the convenience of derivation, we set σ 2 = 0 and α = 4. Based on Equations (10) and (13), we can get the closed-form expression of network EE in the following:

Joint Optimization of FeICIC Parameters and Base-Station Density
Due to the fact that MBSs are usually deployed by network operators, MBS density will change slowly and can be assumed to be constant for analysis simplification. Further simplification of the problem analysis, MBS transmit power in USF and PBS transmit power can also be assumed to be constant without considering power control. Moreover, the PSF ratio can be calculated according to Equation (7). Hence, the network EE is mainly impacted by pico CRE bias, power reduction factor, and PBS density under different UE density, i.e., network load. As MBS density is a constant value, we can first obtain the optimal value of the ratio between PBS density and MBS density λ p,m , denoted as λ * p,m , to maximize the network EE. Then the optimal PBS density λ * p can be calculated by λ * p = λ m λ * p,m . Thus, the joint optimization problem with the object of network EE maximization can be formulated as follows: arg max However, the network EE is nonlinear with λ p,m , B p and ρ, which is difficult to obtain the optimal λ p,m , B p and ρ at the same time with reasonable complexity. Note that the value ranges of λ p,m , B p and ρ are limited, which make it possible to seek out the optimal values of λ p,m , B p and ρ through a linear search algorithm by fixing two of these three variables, respectively. Therefore, we propose a heuristic algorithm to obtain the sub-optimal solution of the joint optimization problem in Equation (15). The proposed heuristic algorithm decomposes the original optimization problem into two sub-problems including FeICIC parameters optimization and PBS density optimization. For FeICIC parameters optimization, we first derive the approximate relation between B p and ρ. Then, we get the sub-optimal values of B p and ρ with given λ p,m by an alternating algorithm. For PBS density optimization, the optimal value of λ p,m can be obtained by a linear search method based on fixed B p and ρ. Finally, we alternately solve two sub-problems to achieve globally optimal values of these variables.

Joint Optimization of Pico CRE Bias and Power Reduction Factor
In order to get the optimal values of B p and ρ for network EE maximization, suppose that λ p,m and λ u are given. Thus, the pico CRE bias and power reduction factor joint optimization problem can be formulated as follow: where B * p and ρ * are the optimal values of pico CRE bias and power reduction factor, respectively. To simplify the solving process, we assume that the overall SINR of PSF MUEs is identical to that of the USF MUEs. Hence, the result of resource allocation will have a direct influence on the network EE. In view of user fairness, the optimal network EE can be achieved when the relationship of association probabilities of PSF MUEs and USF MUEs obey the Equation (17).
Combining with Equation (7), we can get the approximate relation of ρ and B p as: where ρ denotes the approximate near-optimal value of ρ. The relationships between ρ and B p under different λ p are shown in Figure 5. We can find that ρ is a strictly increasing function with respect to B p at a given PBS density. Substituting ρ into Equation (14), we can get the near-optimal value of B p , denoted asB * p by solving the following univariate problem.
Thus,B * p andρ * are obtained by an alternating algorithm, which is shown in Algorithm 1. ∆ 1 and ∆ 2 represent the step lengths of B p and ρ, which are set to be 0.1 and 0.05, respectively. EE * denotes the optimal value of the network EE. The main idea of this altering algorithm is to obtain the near-optimal values of pico CRE bias and power reduction factor by a two-step linear search approach on the basis of the approximate relationship between them. In line 1 of Algorithm 1, the network scenario and relative parameters are initialized. From line 2 to line 10, the optimal pico CRE bias is obtained by a linear search way on the basis of approximate relationship between pico CRE bias and power reduction factor. From line 11 to line 18, the optimal power reduction factor is calculated based on Equation (14). B p = B p + ∆ 1 . 10: end while 11:ρ * = ρ = 0, B p =B * p . 12: while ρ < 1 do 13: Substituting ρ into Equation (14) with given λ u , λ m , λ p,m and B p , the current network energy efficiency EE can be obtained as: EE = EE (ρ) |λ u , λ m , λ p,m , B p . 14: if EE > EE * then 15:ρ * = ρ, EE * = EE . 16: end if 17: ρ = ρ + ∆ 2 . 18: end while

Optimization of PBS Density
Similarly, suppose that B p , ρ and λ u are known. Thus, the PBS density optimization problem can be formulated as follow: Assume that the step length of λ p,m is ∆ 3 , which is set to be 0.03. Thus, the optimal PBS density can be obtained by a linear search algorithm to maximize the network EE, which is described in Algorithm 2. Line 1 of Algorithm 2 indicates the network scenario and parameters initialization. From line 2 to line 8, the optimal PBS density is acquired by a linear search method to maximize the network EE. Substituting λ p,m into Equation (14) with given λ u , λ m , B p and ρ, the current network energy efficiency EE can be obtained as: EE = EE λ p,m |λ u , λ m , B p , ρ. 4: if EE > EE * then 5: λ p,m = λ p,m + ∆ 3 . 8: end while

Joint Optimization of Pico CRE Bias, Power Reduction Factor and PBS Density
The FeICIC parameter optimization sub-problem and the PBS density optimization sub-problem are solved independently by the aforementioned optimization algorithms. Due to the fact that these variables are affected by each other, we further propose a heuristic pico CRE bias, power reduction factor and PBS density joint optimization algorithm to globally optimize network EE based on the joint pico CRE bias and power reduction factor optimization (JBPO) algorithm and PDO algorithm. The detailed procedure of our proposed heuristic algorithm is summarized in Algorithm 3. ε represents the positive tolerance value. N loop is the iteration times of the algorithm. Line 1 of Algorithm 3 signifies the network scenario and parameters initialization. From line 3 to line 4, the JBPO algorithm is executed to obtain the current optimal pico CRE bias and power reduction factor at a given PBS density. From line 5 to line 6, the PDO algorithm is performed to acquire the current optimal PBS density based on the above obtained pico CRE bias and power reduction factor. From line 2 to line 9, the network EE is iteratively optimized until it cannot improve further within an arbitrary value ε. As a result, the optimal pico CRE bias, power reduction factor and PBS density are obtained and the network EE is maximized. Solving the optimization problem in Equation (16), the near-optimal pico CRE biasB * p and the near-optimal power reduction factorρ * at given λ p,m can be obtained based on the JBPO algorithm. 4: Solving the optimization problem in Equation (21), the near-optimal PBS densityλ * p at given B p and ρ can be obtained based on the PDO algorithm. 6: Substituting B p , ρ and λ p,m into Equation (14) with given λ u and λ m , the current network energy efficiency EE can be obtained as: EE = EE B p , ρ, λ p,m |λ u , λ m . 8:

Computational Complexity
The computational complexity of the JBPO algorithm can be calculated as O n B p + n ρ , where n B p and n ρ are the space sizes of the linear search for B p and ρ, respectively. The computational complexity of the PDO algorithm is O n λ p , where n λ p is the space size of the linear search for λ p .
Thus, the computational complexity of the JBPDO algorithm is O n B p + n ρ +n λ p × N loop .
Due to the fact that there does not an exist effective algorithm for solving Equation (15), we compare the computational complexity of our proposed JBPDO algorithm with that of a traversal algorithm. As for solving Equation (15), the optimal values of pico CRE bias, power reduction factor and PBS density can be obtained by a traversal way, i.e., traversal pico CRE bias, power reduction factor and PBS density optimization (TBPDO) algorithm, which refers to traversing all possible values of these three parameters to maximize the objective function in Equation (15). Thus, the computational complexity of TBPDO algorithm will be O n B p × n ρ × n λ p , which signifies that the computational complexity of our proposed JBPDO algorithm is reduced effectively.
As for solving the objective function in Equation (16), the optimal values of pico CRE bias and power reduction factor can also be obtained by a traversal way, i.e., traversal pico CRE bias and power reduction factor optimization (TBPO) algorithm, which refers to traversing all possible values of these two parameters to maximize the objective function in Equation (16). Indeed, TBPO algorithm consists of two nested ergodic sub-processes: (1) traversal pico CRE bias optimization (TBO) algorithm, which is executed by traversing all pico CRE bias to maximize the objective function in Equation (19) under fixed power reduction factor and PBS density; (2) traversal power reduction factor optimization (TPO) algorithm, which is executed by traversing all power reduction factor to maximize the objective function in Equation (20) under fixed pico CRE bias and PBS density. Therefore, the computational complexities of the TBO and TPO algorithms are O n B p and O n ρ , respectively. As a result, the computational complexity of the TBPO algorithm will be O n B p × n ρ , which is obviously higher than that of our proposed JBPO algorithm.

Numerical Results and Analysis
In this section, we not only provide theoretical simulation, but also verify the effectiveness of proposed heuristic algorithms by Monte Carlo simulation. In the theoretical simulation, we considerws a network coverage area within a square region of 1000 m × 1000 m. The deployments of PBS and MBS follow the PPP model and the typical UE is deployed in the origin. The simulation parameters used in this paper are summarized in Table 2. We took the average results from 30 times of network implementations as Monte Carlo simulations results. In each network implementation, the locations of MBSs, PBSs and UEs were modeled as spatial PPP, respectively. Then, the network EEs was calculated for all the different combination of pico CRE bias, power reduction factor, and PBS density values within their value ranges based on wireless channel quality. Finally, the maximum network EE and the optimal pico CRE bias, power reduction factor and PBS density can be obtained by comparing the calculated network EEs under all combinations. Indeed, the results of Monte Carlo simulation, including the maximal network EE, the optimal pico CRE bias, the optimal power reduction factor and the optimal PBS density, were obtained by a traversal way in each network implementation. Meanwhile, the performances of algorithms shown in the following simulation figures were just simulated based on theoretical derived Equation (14). Therefore, Monte Carlo simulation results had the best performance and can be referred to as a baseline for valuing the performances of those algorithms. Table 2. Network scenario parameters.

Parameters Value
Carrier frequency f 2 GHz Total spectrum bandwidth W At first, the network EE performances of the JBPO algorithm were compared with those of the TPO algorithm, TBO algorithm, theoretical simulation and Monte Carlo simulation, as shown in Figures 6 and 7, respectively. As shown in Figure 6, with the increase of PBS density, network EE increased accordingly. With PBS density increase, more UEs were offloaded into the coverage of low power PBSs. Then the distance between the transmitter and the receiver was shorted, resulting in network EE improvement. As shown in Figure 7, with the increase of UE density, the curves of network EE also rose accordingly. As the UE density increased, more UEs can be offloaded into the coverage of low power PBSs by adjusting pico CRE bias and power reduction factor via network EE optimization algorithms. Then network EE can be improved. For the theoretical simulation curve, because the network EE was not optimized and just calculated according to Equation (14) with pico CRE bias and power reduction factor being set to be fixed values 5 dB and 0.25, respectively, so it had the worst performance. The TPO algorithm and the TBO algorithm optimized power reduction factor and CRE bias, respectively. Therefore, the performances of these two algorithms with one parameter optimized were better than that of the theoretical simulation. Our proposed JBPO algorithm can jointly optimize pico CRE bias and power reduction factor together. Therefore, it can further improve the network EE and match the Monte Carlo simulation results very well. Furthermore, the performance of the TPO algorithm was far less than that of the TBO algorithm, which signifies that the influence of pico CRE bias was more than that of the power reduction factor on the network EE, especially in low PBS density. In addition, in Figure 7 the performance gap between the JBPO algorithm and theoretical simulation increases accordingly with the growth of UE density, which further indicates the importance of joint pico CRE bias and power reduction factor optimization for the heavy network load scenario.   Then, the performances of our proposed JBPO algorithm were compared with that of the TBPO algorithm in Figures 8 and 9 from different aspects, respectively. The relationship between network EE and λ p with different λ u are shown in Figure 8. We can see that the performance of our proposed JBPO algorithm was just slightly worse than that of the TBPO algorithm, but the computational complexity of JBPO was much lower than that of the TBPO algorithm. In addition, all curves of network EE increased first and then fell down slightly with PBS density increase, which illustrates that increasing PBS density can improve the network EE significantly within a certain network load. Nonetheless, when the PBS density exceeded a certain limit, further increasing will cause more complex interference and more power consumption, which will result in the network EE deterioration.   The relationship between network EE and λ u with different λ p are depicted in Figure 9. Due to the curves cross with each other under different UE densities, the PBS density should be carefully adjusted according to the network load fluctuation. In addition, for a given PBS density, we can see that the network EE become greater with a higher network load. Referring to Figures 8 and 9 together, although the network EE performance of the TBPO algorithm is just slightly better than that of the JBPO algorithm, the computational complexity of it is O(n B p × n ρ ), which is much higher than that of the JBPO algorithm.
Finally, the network EE performances of the JBPDO algorithm were compared with those of the PDO algorithm, JBPO algorithm with fixed λ p = 10λ m , TBPDO algorithm, and Monte Carlo simulation in Figure 10. In the TBPDO algorithm, pico CRE bias, power reduction factor and PBS density are jointly optimized to maximize network EE by an exhaustive traversal algorithm based on Equation (14). The Monte Carlo simulation results show the maximum network EE within the value range of pico CRE bias, power reduction factor and PBS density at different network load.
As shown in Figure 10, the proposed JBPDO algorithm can obtain better network EE than that of the JBPO algorithm since the JBPDO algorithm further optimizes the PBS density on the basis of the JBPO algorithm. Meanwhile, the accuracy and effectiveness of our proposed JBPDO algorithm are once again verified by Monte Carlo simulation results. In addition, although the network EE of the TBPDO algorithm outperforms our proposed JBPDO algorithm, the computational complexity of TBPDO algorithm is O n B p × n ρ × n λ p , which is far higher than that of the JBPDO algorithm. The convergence of JBPDO algorithm is provided in Figure 11. From Figure 11, we can find that JBPDO algorithm can converge after three iterations. It is proved that the computational complexity of the JBPDO algorithm is much lower than that of the TBPDO algorithm and more suitable for the real-time network.

Conclusions
In this paper, pico CRE bias, PSF power reduction factor and PBS density are jointly optimized to maximize network EE for a two-tier HetNets with FeICIC. First, we derive the closed-form expression of network EE based on stochastic geometry theory. Then, the near-optimal values of pico CRE bias and power reduction factor are obtained by an alternating algorithm based on the equivalence relation between them at a given PBS density deployment. With fixed pico CRE bias and power reduction factor, the PBS density is optimized by a linear search method. Finally, a heuristic algorithm is proposed to optimize the pico CRE bias, power reduction factor, and PBS density jointly for network EE maximization. Extensive simulation results show the accuracy of network EE theoretical deduction and the effectiveness of our proposed low-complexity heuristic algorithm for network EE improvement.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Proof of Lemma 1
Considering a typical UE, it is at a distance r k away from its nearest BS in the k tier, that is to say, there is no BS closer than r k in the k tier. Due to the locations of MBSs and PBSs follow two iid. PPPs, the cumulative distribution function (CDF) of two-point distance r k is where λ k is the BS density of k tier. We can obtain the PDF of the distance by the differential of Equation (A1) over r as follows: Hence, given an arbitrary coefficient value of n, the probability of r p > nr m can be given as: Prob r p > nr m = Prob (no BS closer than nr m |r m ) = Prob r p > nr m , r m = ∞ 0 Prob r p > nr m , r m f r m (r)dr = Further, given two arbitrary coefficient values of n a and n b , the probability of n a r m < r p < n b r m can be expressed as: Prob n a r m < r p < n b r m = Prob r p > n a r m − Prob r p > n b r m = λ m λ m +n 2 a λ p Combing Equation (A4) with Equation (2), the association probabilities of the typical UE can be obtained as Equation (5).

Appendix B. Proof of Lemma 2
According to the Bayes rule, we can obtain the conditional probability as Prob r m |r p > nr m = Prob r m < r,r p > nr m Prob r p > nr m , where Prob r p > nr m is given in Equation (A3). After taking partial derivation with respect to the variable r, the PDF of distance r between a typical UE to its serving BS under the condition of r p > nr m is Further, we can obtain as: f r m |n a r m <r p <n b r m (r) = d dr Prob r m |n a r m < r p < n b r m = d dr Prob r m < r, n a r m < r p < n b r m Prob n a r m < r p < n b r m = 1 Prob n a r m < r p < n b r m d dr Prob r m < r, r p > n a r m −Prob r m < r, r p > n b r m = 2πrλ m Prob n a r m < r p < n b r m exp −πr 2 λ m + n 2 a λ p − exp −πr 2 λ m + n 2 b λ p . (A7) According to Equations (5) and (A7), we can get the PDFs of distance r between a typical UE l to its serving BS in Equation (6).

Appendix C. Proof of Lemma 3
Due to W l and N l can be seen as constant and f l (r) can be obtained based on Equation (6), we can deduce the closed-form expression of the average achievable downlink rate of UE l as: where P l is full transmission power of serving BS of user type l, r l represents the distance between UE l and its serving BS. ρ l is determined as ρ l = ρ when l ∈ {pm} and ρ l is determined as ρ l = ρ when l ∈ {pm, pp}. Otherwise ρ l = ρ l = 1. Considering the complexity of Equation (A8), we can deduce E [log 2 (1 + γ l )] first as Equation (A9).
where B l denotes the association bias of the typical UE l to select its serving BS. When l = pp, B l = B p ; otherwise B l = 1. In order to simplify the analysis, we can define variable β = B k,l /B l , k ∈ {m, p}, where B m,l =B m / ρ l ρ −1 l and B p,l = 1, if l = up B p , otherwise .
Substituting Equations (A10) and (A11) into Equation (A9), and further substituting Equation (A9) into Equation (A8), we can obtain the closed-form expression of the average achievable downlink rate of a typical UE l as Equation (9).