Planning capacity for 5G and beyond wireless networks by discrete fireworks algorithm with ensemble of local search methods

In densely populated urban centers, planning optimized capacity for the fifth-generation (5G) and beyond wireless networks is a challenging task. In this paper, we propose a mathematical framework for the planning capacity of a 5G and beyond wireless networks. We considered a single-hop wireless network consists of base stations (BSs), relay stations (RSs), and user equipment (UEs). Wireless network planning (WNP) should decide the placement of BSs and RSs to the candidate sites and decide the possible connections among them and their further connections to UEs. The objective of the planning is to minimize the hardware and operational cost while planning capacity of a 5G and beyond wireless networks. The formulated WNP is an integer programming problem. Finding an optimal solution by using exhaustive search is not practical due to the demand for high computing resources. As a practical approach, a new population-based meta-heuristic algorithm is proposed to find a high-quality solution. The proposed discrete fireworks algorithm (DFWA) uses an ensemble of local search methods: insert, swap, and interchange. The performance of the proposed DFWA is compared against the low-complexity biogeography-based optimization (LC-BBO), the discrete artificial bee colony (DABC), and the genetic algorithm (GA). Simulation results and statistical tests demonstrate that the proposed algorithm can comparatively find good-quality solutions with moderate computing resources.

MIMO, enable device-to-device (D2D) communication, reduce the protocol overhead, implement heterogeneous cell architecture, allow network function virtualization, install content caching, and improve the range of spectrum for a typical 5G and beyond wireless network [4].
The seamless planning of the 5G and beyond wireless networks is a one-time activity before a system is initially installed. Therefore, we aim to incorporate only necessary high-level network details in our proposed 5G and beyond wireless network model to keep its design simple. However, heterogeneity is one of the main characteristics of the 5G and beyond wireless networks, which means various wireless network technologies work in tandem with robust software and hardware implementation. Although, several other choices are available to define the link-rate between two wirelessly communicating nodes [5][6][7][8][9][10][11]; however, for the simplicity and to incorporate only necessary high-level network details, we use path loss as a criterion of variation for data rates between a wireless link of two communicating nodes. A rationale behind using the path loss as a data rate variation criterion is to cater to the heterogeneous and multi-RAT characteristics of the 5G and beyond wireless networks. Typically, the goal of WNP is to minimize the overall system operating cost. Such cost includes the BSs costs, RSs costs, installation, and operational cost of a network. Here, path loss is included as an operational cost among each communicating link. Also, we incorporate the practical network constraints such as maximum load on BSs and RSs.
The model presented in this paper allows single-hop communication only because this model is specifically used to plan the capacity of a network in the densely populated urban/commercial centers. The proposed model can be used only for planning 5G and beyond wireless networks from scratch. This work aims to devise a methodology to plan simultaneous BSs and RSs locations for the 5G and beyond wireless networks.

Contributions
In this paper, we propose a relatively new swarm intelligence-based algorithm to solve the proposed WNP plan. Following are the main contributions of this paper: 1. Two equivalent integer programming problems are formulated for the 5G and beyond wireless network planning: • 5G and beyond wireless network planning as a binary search space.
• 5G and beyond wireless network planning as an integer search space.
2. The proposed planning problem is solved by a discrete fireworks algorithm (DFWA) using ensemble of local search methods: "swap," "interchange," and "insert." • A repair algorithm is proposed to repair the illegal candidate solutions.
• A T test is conducted to statistically analyse the performance of the proposed algorithm.
In the rest of the paper, the existing work related to planning is discussed in Section 2, and the system model and mathematical framework are presented in Section 3. The proposed Firework-based algorithm is discussed in Section 4, and its computational complexity is presented in Section 5. The simulation results are discussed in Section 6, and finally, the paper is concluded in Section 7.

Related work
Conventionally, there are two types of WNP: by coverage and by capacity. Here, multiple RSs can be deployed between user equipment (UE) and a BS while planning the coverage of the 5G and beyond wireless network. In contrast, only a single RS can be deployed between a UE and a BS to enhance the capacity of a network. In literature, various models for the coverage of the 5G and beyond wireless network is presented [6,7,12]. In this paper, we discuss capacity planning for wireless networks.
Deploying an RS within the serving area of a BS can increase network throughput, however, this deployment raises the overall hardware cost. The authors present a deployment algorithm for an IEEE 802.16j network. In [9], a three-phase RS deployment algorithm is proposed. The first phase aims to construct several promising zones where an RS can be deployed. An RS deployment in each zone will improve the transmission rate from a subscriber station (SS) to a BS. In the second phase, larger zones are constructed by combining several smaller zones to reduce the number of deployed RSs. When all the RSs are deployed in promising zones, the results show that the transmission delay and the hardware cost is reduced by using the proposed algorithm. In [8], RS location planning is formulated for capacity gains in the IEEE 802.16j network transparent mode. The RS locations are determined from a given set of candidate RS sites, and the optimization problem is solved by the interference aware algorithm. An integer programming problem is addressed to determine BS and RS locations that will enhance network capacity at a minimal cost. A similar problem for the simultaneous BS/RS problem is formulated to determine BS and RS locations to minimize the operational and hardware cost [10,11].
In [5], a 4G/5G heterogeneous network (HetNet) is proposed for small cells (SCs) with additional feature of fault-tolerance. SCs are intended to increase network capacity and to extend network coverage (i.e., multi-hop) for better spectrum efficiency. An expanded approach is adopted to avoid nonlinearity in the mathematical model of a network. In [1], two infrastructure-aware planning strategies are proposed for SCs and fiber back-haul. These strategies include joint design (JD) for the cost minimization of SCs with fiber backhaul and traditional design (TD) based on a two-step optimization approach. The relative performance of these two strategies is compared.
In [6], a 5G network model is presented to deploy infrastructure with machine type communication. The goal of this research is to plan 5G infrastructure with the best coverage under various constraints. The formulated problem is multi-objective integer programming problem, and heuristic algorithms are presented to solve this challenging optimization problem. In contrast, a mixed integer linear programming (MILP) problem for the 5G network is formulated in [13]. The goal of the optimization problem is to maximize operating profit of the proposed 5G network. The problem is solved using the CPLEX software package. Note that CPLEX software package is presented by the IBM and this software package implements optimizers based on the simplex algorithms [14].
The multi-tenancy model is used in the fifth generation of mobile networks in which diverse operators share the same wireless infrastructure. To achieve such goals [15], small cells (SCs) offer network providers more flexible, scalable, and cost-effective solutions compared to the macro-cell deployments. The authors proposed a framework for cell planning in multi-tenant SC networks that is evaluated considering the deployment of a new tenant, where different sets of planning specifications are tested. Similarly, in [16], a framework for automated cell planning in multi-tenant SC networks is investigated. By taking advantage of the available network data, a set of detailed planning specifications over time and space domains are generated to meet the improved capacity by each tenant. Then, the network infrastructure and configuration are updated by using an algorithm that considers different actions such as adding/removing channels and adding or relocating SC.
In [17], nomadic nodes (NNs) enabled 5G network deployment is proposed to provide demand-driven service, which not only increases the network capacity and extends the cell coverage, but also reduces network energy consumption. The 5G network benefits from location information, which is used in wireless network design and optimization. Several challenges can be addressed by using location information such as an increase in data traffic, accommodate an additional number of devices, robustness in mission-critical services, reduction in total energy consumption, and improvement in latency [18]. Reliability and latency are some of the sensitive issues that need to be considered during node deployment in the 5G networks. To ensure network reliability and latency, key resources such as computing, network, and storage need to be brought at the edge of the 5G network. In [19], such a key scenario is taken into consideration to achieve robustness and cost-effectiveness while planning the capacity of the 5G network.
The objective of the 5G and beyond wireless network planning is to minimize infrastructure (BSs and RSs) and operational costs. The proposed network plan is to enhance the capacity of the network in densely populated urban centers. This work has two distinctions when compared against exiting work. First, two types of mathematical problem formulations are presented with the goal of reducing the number of constraints checks during implementation, which may help to further explore mathematical aspects of the problem. Second, a new fireworks algorithm with ensemble of local search methods is presented to solve the planning problem.

System model and mathematical framework
The model of the proposed 5G and beyond wireless network with the aim of capacity planning is shown in (Fig. 2). The proposed network is comprised of three types of nodes: user equipment (UE), relay stations (RSs), and base stations (BSs). In this section, the system model of the proposed 5G and beyond wireless network is briefly discussed.

System model
In the proposed 5G and beyond wireless network, a UE can be a subscriber or a group of subscribers with certain data traffic demand. A UE can be connected to a BS directly or communicated to a BS via an RS subject to fulfill its traffic demand. UEs may or may not be battery powered. A UE cannot communicate through more than one RS or more than one BS or both. An RS is used to relay data from a BS to UE subject to fulfill the given constraints. Here, we only consider data traffic in the downlink direction in the proposed network model. RS is not directly connected to the Internet. However, a BS is directly connected to the Internet and has reasonable computing power. An RS has less computing power as compared to a BS. Moreover, one or more RSs can be connected to a BS.

Problem formulation
The objective of the 5G and beyond wireless network planning is to minimize the overall cost (hardware, operational) of network functioning in the presence of users' traffic demand. Symbols used in this paper are described in ,t denote path loss and upper-bound (i.e., channel capacity) on the possible information flow-rate associated with the links (b, r), (r, t) and (b, t) respectively. Note that u T t denotes the traffic demand of a user equipment t.

Cost function
Low hardware and operational costs are part of the objective in the proposed planning to develop a 5G and beyond wireless network. Total hardware expenses include the cost of deployed BSs and deployed RSs at their respective selected sites. UEs can communicate to a deployed BS directly or indirectly via a deployed RS, and a deployed RS must communicate to a deployed BS. Communication among 5G and beyond wireless network nodes (UE, BS, RS) occurs at a lower cost when nodes are in close proximity. Cost function consists of two parts, i.e., hardware cost and operational cost. Total cost of the proposed planning problem is the weighted sum of the hardware and operational costs. In communication among distant nodes, more power is consumed, the quality of the communication may be degraded, and the wireless links might be impaired due to path loss. Therefore, we incorporate path loss as an operational cost for each communicating link (i.e., user equipment (UE)-BS, UE-RS, RS-BS). The cost function and constraints for the 5G and beyond wireless network planning are as follows: such that:

Topology constraints
x RT r,t ≤ y R r , ∀r ∈ R, ∀t ∈ T.

Flow constraints
The flow f BT b,t is a function of binary decision variable and is defined as follows: , the flow value f BT b,t must be within the maximum capacity m BT b,t as follows: The flow f RT r,t is a function of binary decision variable and is defined as follows: For each existing RS-UE connection (i.e., x RT r,t =1), the flow value f RT r,t must be within the maximum capacity m RT r,t as follows: The flow f BR b,r is a function of binary decision variables x BR b,r , x RT r,t and is defined as follows: , the flow value f BR b,r must be within the maximum capacity m BR b,r as follows:

Load constraints
In the cost function given in Eq. (1), the first term represents the used hardware and installation costs of a BS and an RS, respectively. As communication between nearby nodes (UE, BS, RS) is more promising with respect to the cost and clarity than communication between distant nodes. A second term is used in Eq. (1) to incorporate the notion of nearby communication. The second term in (1) presents the path loss of each established communication link between communicating nodes. Here, W 1 and W 2 are weight parameters for the first and second terms, respectively, as shown in Eq. (1). The topology constraint (2) ensures that each RS, if deployed, is connected to one BS only, and constraints (3) -(4) define that every UE and RS can be connected only to a deployed BS. Constraint (5) confirms that every UE is connected to a deployed RS. Constraint (6) ensures that each UE is either connected to a BS or an RS, but not both. Constraints (7) -(9) ensure that the flow value on links (b, r), (b, t), (r, t) is within maximum capacity. Constraints (10) and (11) confirm that the load on each deployed BS and RS does not exceed the corresponding maximum load. Finally, (12) guarantees that every UE has enough flow through either a BS or an RS.

Redefining decision variables
In the wireless network planning problem R, B, and T denote sets RS sites, BS sites, and UE, respectively. The proposed network planning problem has an integer domain, and we can define a candidate solution as a vector of non-negative integers as follows: where |T| and |R| are the cardinalities of sets T and R. Note that in X, T t represents the ith UE connected to some BS b or RS r , and R r represents the rth RS that can be connected to some BS b. Here, R r is zero when the rth RS is not deployed. In X, each UE can be connected to anyone BS or RS. Therefore, in the candidate solution, we represent BS, followed by RS in consecutive orders of positive integers. Suppose, we have |B| BSs and |R| RSs sites in a candidate solution X; then, in the candidate solution X representations for BSs sites are 1, 2, . . . , |B| and representations for RS sites are 1 + |B|, 2 + |B|, 3 + |B|, . . . , |R| + |B|.
In X, integer variable T t , t = 1, 2, . . . , |T| takes a value in set {1, 2, . . . , |B| + |R|}, where element 1, 2, 3, . . . , |B| represents the BS sites, and |B| + 1, |B| + 2, |B| + 3, . . . , |B| + |R| represents the RS sites. Integer variable T t indicates UE t is connected to a BS or a RS site. T t =i if i ≤ |B| indicates that UE t is connected to a base station site i and a base station is installed at the base station site i. T t =i if i>|B| indicates that UE t is connected to relay station site i − |B| and the relay station is deployed at site i − |B|. Integer variable R r takes a value in 0, 1, 2, . . . , |B|. Here, R r =0 indicates that a relay is not installed (deployed) at relay site r. Moreover, R r =0 also indicates the RS is not connected to any UE or any BS; hence RS is not deployed. R r =b, if 1 ≤ b ≤ |B| indicates that a relay is installed (deployed) at relay site r and that the relay site is connected to a base station at site b. The same also indicates that BS site b should have a base station installed (or deployed). Both T t and R r variables decide whether there is a BS in BS site b. The parameters of the WNP problem are the same before and after encoding, except a small adjustment for the experimentation. Now, we reformulate the WNP problem using variable X; the indices are used: t = 1, 2, . . . , |T|, r = 1, 2, . . . , |R|, and b = 1, 2, . . . , |B|.

Cost function
We reformulate the plan for the WNP problem using the encoded vector X. The cost function of the WNP is used to minimize the overall cost of the network. Reformulations of the cost function and constraints for the WNP are given as follows: Let α B b represents whether a BS is in X: Let β R r represent whether a RS is in X: such that:

Topology constraints
Making use of the topology constraints (2) -(6) from the formulation and to reduce the number of variables, we define a vector of non-negative integers X = (T 1 , T 2 , T 3 , ...T |T| , R 1 , R 2 , R 3 , ..., R |R| ), where |T| and |R| are the cardinalities of sets T and R. Constraints (2) -(6) are implicitly enforced in the vector X.

Flow constraints
The flow f BT b,t is a function of decision variable X and is defined as follows: For each connection between a BS and a UE, the flow value f BT b,t must be within the maximum capacity m BT b,t as follows: The flow f RT r,t is a function of decision variable X and is defined as follows: The flow f BR b,r is a function of decision variables X, and is defined as follows: , if X j = b, j = r + |T|, ∀b = 1, 2, ...|B|, ∀r = 1, 2, ...|R|

0, Otherwise
For each existing connection between a BS and an RS, the flow value f BR b,r must be within the maximum capacity m BR b,r as follows:

Load constraints
Constraints (18) and (19) ensure that the load on each BS and RS does not exceed the maximum load and (20) guarantees that every UE has enough flow, either through a BS or an RS, but not both. Redefining the decision variables in the Section 3.3 enforce some constraints of the planning problem implicitly, which reduces the number of constraints checks during implementation. Therefore, second formulation has reduced number of constraints checks.
The proposed WNP is a combinatorial integer space optimization problem, and to the best of author's knowledge, no known polynomial-time algorithm exists to solve such problems. The exhaustive search requires high computing costs to find an optimal solution for such problems. A practical approach is to use approximate algorithms such as evolutionary algorithms (EAs) to solve challenging problems in moderate computing resources. Without any guarantee of an optimal solution, EAs can provide a high-quality solution using moderate computing resources. Therefore, we propose a new swarm intelligence-based EA, i.e., discrete fireworks algorithm (DFWA), with an ensemble of local search to solve the formulated wireless network problem.

Swarm intelligence-based evolutionary algorithms
Like the biogeography-based optimization (BBO) [14,[20][21][22][23] and the discrete artificial bee colony (ABC) [23] algorithms, discrete Firework algorithm (DFWA) is a swarm intelligence-based algorithm [24,25]. Individuals in these population-based metaheuristics act as a simple agent and behave collectively in a decentralized, self-organized manner. Individual agents in a typical swarm intelligence-based algorithm can communicate either directly or indirectly with each other by acting in its local environment [26]. An individual agent of a swarm follows straightforward rules; however, interactions between such agents can be complex, causing global behavior that is far beyond the capability of the individual agents.

DFWA with an ensemble of dynamic local search methods
Selecting a local search (LS) method for any evolutionary algorithm is a tedious task, especially when uncertainty in the performance of the LS operators is expected [25]. In the absence of experimental evidence and in the presence of variations in the performance of LS methods, a dynamic ensemble of local search methods to avoid manual selection of LS methods that might lead to unsatisfactory performance. Therefore, inspiring from [25], we present a DFWA with an ensemble of dynamic LS methods (DFWA-with-Dy-3-LS) and LS methods include in the ensemble are insert, interchange, and swap. These LS methods are defined and discussed in [14,25]. Like the DFWA [24], the DFWA-with-Dy-3-LS algorithm has four parts: an explosion operator, a mutation operator, a repair mechanism, and a selection strategy.

Explosion operator
In the DFWA-with-Dy-3-LS algorithm, the cost value and the control parameters determine the criteria of the explosion operator. The explosion operator uses LS methods with two control parameters: explosion strength and explosion radius. The DFWA-with-Dy-3-LS explosion operator determines the number of sparks and the radius of those sparks in proportion to the cost value of fireworks.
Explosion strength: In the DFWA-with-Dy-3-LS, the explosion strength determines the number of sparks that are generated by the explosion of a firework. The cost of a firework and user-defined control parameters determine the number of sparks that are generated by a firework. Like the DFWA, the DFWA-with-Dy-3-LS is designed in such a way that a firework with a lower cost (good firework) generates more sparks around that firework. A firework with a higher cost (bad firework) should generate fewer sparks around that firework. The rationale behind generating more sparks around a good firework is to exploit the low cost of the good firework, and a thorough search is conducted to find a better solution around the good firework. However, a bad firework should generate sparks that are fewer in number and sparse in distribution to avoid unnecessary computing. The sparks generated from the bad fireworks are used to explore the search space and to prevent the algorithm from being trapped in a local minimum. The DFWA-with-Dy-3-LS computes the number of sparks,s i , for the ith firework, where i = 1, 2, . . . , N.
where s i is the number of sparks for the ith firework (for each of i = 1, 2, ..., N), y max is the maximum cost of the N fireworks in the current algorithm generation, f (X i ) represents the cost of the ith firework, M e is a constant that controls the total number of sparks generated by N fireworks, and is a small constant used to avoid a division by zero in (21).
Dynamically selecting an LS method from an ensemble of LS methods: In the DFWAwith-Dy-3-LS, we propose a new criterion to dynamically select an LS method from an ensemble of LS methods during DFWA operation. We denote three LS methods as a set of integers ℘. For example, the set ℘ = {1, 2, 3} represents an ensemble of LS methods in which elements of the set ℘ represent the three LS methods Op1=1, Op2=2, and Op3=3, as insert, interchange, and swap respectively. Initially, the DFWA-with-Dy-3-LS algorithm randomly assigns an LS method from the set ℘ to each of the N fireworks. In the first generation of the DFWA-with-Dy-3-LS, each firework generates, s i , sparks using randomly assigned LS methods, where i = 1, 2, . . . , N, and the generated sparks are evaluated using the cost function (14). In the second generation of the DFWA-with-Dy-3-LS, if there is no improvement observed in the cost value of the sparks generated from the ith firework for each of i = 1, 2, . . . , N fireworks, then, the currently assigned LS method is replaced for the ith firework in the next algorithm generation. The ith firework is randomly selected from the remaining list of LS methods. For example, if ℘ 1 , ℘ 2 and ℘ 3 are elements of set ℘ of the three LS methods, then ℘ 1 is associated with the X 1 firework. In any DFWAwith-Dy-3-LS generation, if there is no improvement observed in the sparks generated by X 1 firework, then ℘ 1 is replaced by randomly selecting an LS method from ℘ 2 and ℘ 3 .
If the population of fireworks is N =10 in the DFWA-with-Dy-3-LS algorithm, then each component of the integer vector ℘ = (3, 3, 1, 2, 2, 3, 3, 1, 1, 1) represents the LS methods associated with the corresponding fireworks. For example, in the population of N fireworks, the first firework is associated with the third LS method, while the last firework is associated with the first LS method.
Explosion radius: The explosion radius is an integer value used to determine the number of times a local search (LS) operator is applied to perturb one or more components of the ith firework. The cost of the ith firework, for each of i = 1, 2, . . . , N fireworks, and control parameters are used to determine the number of times an LS method is applied on that firework [25]. In the DFWA-with-Dy-3-LS, a firework with a lower cost should generate sparks with a smaller radius around that firework, and a firework with a higher cost function value should generate sparks with a larger radius around that firework. The rationale behind generating sparks with a smaller radius is to exploit the good firework, and a thorough search is conducted to find a better solution around the good firework. However, sparks generated from the bad fireworks with larger radius explore the search space and prevent the algorithm from being trapped in a local minimum. The DFWA-with-Dy-3-LS computes the explosion radius, s i , for the ith firework as: where H i is the explosion radius associated with the ith firework (for each of i = 1, 2, ..., N), y min is the minimum cost of the N fireworks in the current algorithm generation, f (X i ) represents the cost of the ith firework,â is a constant that controls the maximum number of times an LS operator is imposed on X i , and is a small constant used to avoid a division by zero in (22).

Mutation operator
In the DFWA-with-Dy-3-LS, a modified mutation operator that uses the random integer function randi for the mutation is adopted. A set of fireworks Z for mutation are selected from N fireworks to set up sparks with the mutation operator, where |Z| < N and |Z| is the cardinality of set Z. One or more components of a firework X i ∈ Z are probabilistically selected with the user-determined probability (2020) 2020:185 Page 14 of 24 mutateProb and replaced with the new component. The DFWA-with-Dy-3-LS generates one spark for each mutation spark X i ∈ Z using the mutation operator as follows: where X i q is the component of a newly generated spark to replace X i q in the current algorithm generation and X min q and X max q are lower and upper bounds of the search space in dimension q.

Repair mechanism
The candidate solutions may become infeasible due to violation of rectangular or nonrectangular constraints during the operation of evolutionary algorithms (EAs). These infeasible solutions are useless for further evolution in any EA. Similarly, randomly generated fireworks and sparks may fall in the infeasible space after executing the DFWAwith-Dy-3-LS operations. Sparks in the infeasible space are considered illegal and are useless for further algorithm operation. Therefore, illegal sparks need to be returned to the feasible space. A candidate solution, as defined in (13) of the WNP, is illegal if it falls outside the feasible space. Therefore, we proposed a repair algorithm for the WNP that has non-rectangular constraints.
Repair algorithm: The pseudo-code for the repair algorithm is presented in Algorithm 1. In the repair algorithm, a candidate solution, X in (13), that is generated randomly or is evolved by an evolutionary procedure, may violate one or more constraints of the optimization problem. Violation of any constraint makes the candidate solution illegal. In the DFWA-with-Dy-3-LS, the legality of each randomly generated firework is verified using the repair algorithm. Further, in each DFWA-with-Dy-3-LS generation, the legality of all the sparks generated using the explosion operation, and mutation operation is verified using the repair algorithm.
Initially, we set the parameters c B b , c R r , l BR b,r , l BT b,t , l RT r,t , m BR b,r , m BT b,t , m RT r,t , u T t , C 1 , C 2 , W 1 , and W 2 according to user-defined criteria. For each candidate solution, X in (13), of the population, we compute f BR b,r , f BT b,t , f RT r,t and check that the flows associated with each link do not exceed the maximum link capacity (i.e., m BT b,t , m RT r,t , m BR b,r ). We consider a link to be illegal if it violates the maximum link capacity, i.e., We disconnect UE and RS from illegal links and reconnect them to legal corresponding links, i.e., UE-BS, UE-RS, and RS-BS: Next, we compute the load for each BS and RS in X. We disconnect UE from overloaded BS and RS and reconnect them to BS and RS subject to the load constraints C 1 , C 2 and legal links f BT b,t < m BT b,t , f RT r,t < m RT r,t . Our proposed iterative process to repair the candidate solution X in (13) does not guarantee that each illegal solution has become a legal solution using limited computing. Here, limited computing means that we are not exhaustively checking each UE-BS, UE-RS, and RS-BS connection to verify its legality. If a candidate solution is not repairable within limited computing, we randomly generate a candidate solution X and check its legality. This process of randomly generating a solution continues until we generate a legal candidate solution.  Compute link flows among TP-RS, TP-BS, BS-RS links. 3 Disconnect an infeasible BS-RS link try to establish a feasible link between a BS-RS. 4 Disconnect an infeasible BS-TP link and try to establish a feasible link between a BS-TP/RS-TP. 5 Disconnect an infeasible RS-TP link and try to establish a feasible links between a BS-TP/RS-TP. 6 Update the candidate solution X. 7 end 8 while links in X is still infeasible do 9 Randomly generate a candidate solution X. 10 Repeat steps 1 to 7. 12 Compute load on each BS and RS in X. 13 for load on each BS/RS in X do 14 if load on a BS/RS is infeasible then 15 Disconnect TPs from an overloaded BS/RSs and reconnect to an underloaded BS/RS subject to maximum link capacity and maximum loads on BSs and RSs. 16 Update the candidate solution X. 17 end 18 end 19 while load on BSs and RSs is still infeasible do 20 Randomly generate a candidate solution X. 21 Repeat steps 1 to 18.

Selection strategy
Each generation of the DFWA-with-Dy-3-LS produces a number of candidate solutions greater than the N fireworks population. Therefore, after applying all the DFWA-with-Dy-3-LS operators, a new N fireworks population needs to be selected from the current candidate solutions. The DFWA-with-Dy-3-LS adopts the same elitism-random selection strategy as that adopted in the enhanced fireworks algorithm (EFWA) [24]. In the DFWA-with-Dy-3-LS, first, the solution with the least cost value is selected, then (N − 1) candidate solutions are randomly selected from the remaining candidate solutions for the next algorithm generation. The pseudo-code for the DFWA-with-Dy-3-LS is given in Algorithm 2. : The best solutions X found so far. Initialization: Randomly generate population of N fireworks; X i where i = 1, 2, ..., N; Randomly generate an integer vector ℘ of size N; Empty set S 1 for candidate solutions.
1 Check the legality of the N fireworks using the repair algorithm, as algorithm 1 pseudo code, and evaluate using the cost function in (14).
2 while stopping criteria not satisfied do 3 for each candidate solution X in population of size N do 4 Calculate s i using formula in equation (21) 5 for j ← 1 tos i do 6 Get a copy of firework to generate sparks X j = X i .

7
Select the ith local search method from the vector ℘ .

8
Calculate H i H i is calculated using equation (22) 9 for k ← 1 toH i do 10 I mpose ith local search method (from ℘ ) on X j . Randomly select a set Z of fireworks for mutation from the population of the N fireworks, where N>|Z|. 17 for each firework X i ∈ Z do 18 Get a copy of firework to generate sparks X j = X i .

19
for q ← 1 tod do 20 if rand < mutateProb then 21 Update X q j using equation (23). Check the legality of the S 1 solutions using the repair algorithm, as algorithm 1 pseudo code, and evaluate using the cost function in (14). 27 Select the best solution and the (N − 1) solutions to make a new population of N fireworks for next algorithm iteration.

28
return the best solution found so far.

Computational complexity
Typically, the computational complexity of population-based evolutionary algorithms (e.g., genetic algorithm (GA)) is analyzed in terms of the number of cost function evaluations [23,27], as in our case (14). However, the computational complexity is highly dependent on the coding efficiency. In our experimental algorithms, the low complexity biogeography-based optimization (LC-BBO) [27], and the GA, the cost function evaluations are equal to GN, where G is the total number of algorithm iterations and N is the population size. As in the LC-BBO algorithm and the GA, the cost function is usually evaluated for a candidate solution at least once in an algorithm generation. However, the cost function evaluation may be made more than once for a candidate solution; the cost function evaluation is made more than once in discrete artificial bee colony (DABC) [23] and DFWA-with-Dy-3-LS.

DFWA-with-Dy-3-LS
Similar to the DABC algorithm [23], the DFWA-with-Dy-3-LS run cost function evaluations more than once in one algorithm generation. Initially, in the DFWA-with-Dy-3-LS, the population of N fireworks is evaluated using the cost function. Then, for each firework i = 1, 2, . . . , N, the number of sparks, s i , generated using explosion operations is evaluated using the cost function (14). We denote M e as the total number of sparks (or candidate solutions) generated during the explosion operation as follows: The DFWA-with-Dy-3-LS selects a set of Z fireworks for mutation from the N fireworks to set up sparks by the mutation explosion, where |Z| < N and |Z| is the cardinality of the set Z. For each firework for mutation, X i ∈ Z, the total number of |Z| sparks are generated and are evaluated using the cost function. Initially, the population of N fireworks is evaluated using the cost function. Then, in each DFWA-with-Dy-3-LS generation, the total number M e and |Z| sparks are evaluated using the cost function. Here, the total number of cost function evaluations in one algorithm generation would be [14]: 6 Results and discussions

Simulation setup
The simulation was run using MATLAB 2017 on i7 series processor with 8 GB of RAM. For each problem instance, the experiment requires inputs such as number of BSs, number of RSs, and number of UEs. The experiment was conducted with eight different problem instances. Problem specific parameters such as the number of BSs, RSs, and UEs are shown in Table 2, and algorithm-specific parameters are shown in Table 3. The UE demand is a real vector which is randomly generated in the interval [0.01 4.0]. The cost for each installed BS and RS is set as BS = 25 and RS = 5. The real matrices representing the path loss for each link l BR b,r , l BT b,t and l RT r,t were randomly generated. The maximum link rates m BR b,r , m BT b,t and m RT r,t for each existing link (i.e., x BR b,r = 1, x BT b,t = 1, and x RT r,t = 1) is defined in Table 4. In this paper, we use discrete artificial bee colony (DABC) algorithm from [23] and BBO's low-complexity version (i.e., LC-BBO) is taken from [22].

Performance
The average cost of the WNP is plotted for the DFWA-with-Dy-3-LS, the LC-BBO algorithm, the DABC algorithm, and GA in Fig. 3. The DFWA-with-Dy-3-LS algorithm outperformed the LC-BBO, the discrete ABC algorithms, and the GA in terms of average cost. On the other hand, the LC-BBO algorithm performs better in terms of average cost against the discrete ABC algorithm and the GA, as shown in the Fig. 3. However, the GA is the worst performer in terms of average cost against all the experimented algorithms, as shown in the Fig. 3. The standard deviation (Std.) for the genetic algorithms (GA) is smaller than the standard deviation of the LC-BBO, the discrete ABC, and the DFWA-with-Dy-3-LS. The LC-BBO algorithm has a smaller Std. than the DFWA-with-Dy-3-LS algorithm. However, DFWA-with-Dy-3-LS has a higher Std. than all the experimented algorithms (see Fig. 4). The DFWA-with-Dy-3-LS algorithm consumed a higher average Matlab CPU time than the LC-BBO, discrete ABC, and GA. The GA consumed the lowest average CPU time as compared to all the experimented algorithms (see Fig. 5).

Performance significance of the DFWA-with-Dy-3-LS
A T test shows a significant difference in the performance of the DFWA-with-Dy-3-LS algorithm, LC-BBO, discrete ABC, and genetic algorithm (GA). The null hypothesis H 0 states that both algorithms produce the same average cost. Also, we performed the T test of an alternative hypothesis H 1 , which states that the DFWA-with-Dy-3-LS produces a lower average cost. Table 5 shows the p values of the T test for each problem instance against each compared algorithm. The p values can be compared against the generally acceptable level of significance α = 0.05 to decide whether hypothesis H 1 is accepted. If  the average cost by the DFWA-with-Dy-3-LS is lower than any compared algorithm and p ≤ α, then we conclude that there is a statistically significant difference between the DFWA-with-Dy-3-LS and the other experimented algorithms. Otherwise, we conclude that the observed difference is not statistically significant.
The DFWA-with-Dy-3-LS showed a lower average cost when compared to the other experimented algorithms, and the p value was also lower than 0.05. Therefore, the performance of the DFWA-with-Dy-3-LS was significantly better than the performance of the LC-BBO algorithm, discrete ABC algorithm, and GA.

Performance analysis
A boxplot in Fig. 6 shows the comparative performance of the algorithms tested in the formulated WNP problem. Although variability and consistency are observed among the experimented algorithms, as shown in Fig. 6. For example, the DABC algorithm and the GA have low variability in the depicted boxplots. The better performance of the DFWA-with-Dy-3-LS is observed in all the experiments, as the DFWA-with-Dy-3-LS achieves better average cost value as compared to the LC-BBO algorithm, discrete ABC algorithm, and genetic algorithm. Using boxplots (in Fig. 6), we can see the outliers (the red "+") in the data set. For example, the DFWA-with-Dy-3-LS, discrete ABC, and genetic algorithm data are symmetric, as shown in the Fig. 6. However, left and right skewness is observed for the LC-BBO and the DFWA-with-Dy-3-LS algorithms, respectively.

Conclusion and future work
In this paper, we propose a wireless network planning (WNP) for 5G and beyond networks as an integer programming problem with a single-hop configuration. This network planning problem consists of three types of nodes: base stations (BSs), relay stations (RSs), and user equipment (UE). A UE can communicate with a BS directly or via an RS. We use path loss as a criterion of variation for data rates between a wireless link of two communicating nodes to cater to the heterogeneous and multi-RAT characteristics of the 5G and beyond wireless network. The objective of this WNP was to minimize the overall hardware and operating cost of the network. Finding an optimal solution using exhaustive search was impractical due to the high computing demand. To the best of our knowledge, no known polynomial-time algorithm exists to solve such computationally challenging problems. Due to performance differences in LS methods of the DFWA, randomly selecting an LS method could result in an inefficient choice; therefore, we proposed a DFWA with an ensemble of three dynamic LS methods (DFWA-with-Dy-3-LS). The DFWA-with-Dy-3-LS considers performance as a metric to select LS method from ensemble during the operation of the algorithm. We experimentally compared the performance of the proposed DFWA-with-Dy-3-LS algorithm against the low-complexity BBO, discrete ABC, and genetic algorithms. The DFWA-with-Dy-3-LS algorithm outperforms the other algorithms in terms of the average cost of the network and is significantly better than the other experimented algorithms. Statistical analysis showed that the DFWA-with-Dy-3-LS algorithm performed significantly better than the low-complexity BBO, discrete ABC, and genetic algorithms. Our experimental results also demonstrate that low average cost and low CPU time can be used to select an appropriate algorithm for the planning of 5G and beyond wireless networks.
We are planning to extend this work in different dimensions in the future. First, we would like to prove/disprove NP-completeness of this mathematical problem. Second, we would like to extend this work to implement some exact algorithms and their performance would be compared with the current approximate algorithms. Third and final, we would like to not only experiment DFWA with various local search methods individually but would also like to expand the ensemble size.