Backhaul-Aware Optimization of UAV Base Station Location and Bandwidth Allocation for Profit Maximization

Unmanned Aerial Vehicle Base Stations (UAV-BSs) are envisioned to be an integral component of the next generation Wireless Communications Networks (WCNs) with a potential to create opportunities for enhancing the capacity of the network by dynamically moving the supply towards the demand while facilitating the services that cannot be provided via other means efficiently. A significant drawback of the state-of-the-art have been designing a WCN in which the service-oriented performance measures (e.g., throughput) are optimized without considering different relevant decisions such as determining the location and allocating the resources, jointly. In this study, we address the UAV-BS location and bandwidth allocation problems together to optimize the total network profit. In particular, a Mixed-Integer Non-Linear Programming (MINLP) formulation is developed, in which the location of a single UAV-BS and bandwidth allocations to users are jointly determined. The objective is to maximize the total profit without exceeding the backhaul and access capacities. The profit gained from a specific user is assumed to be a piecewise-linear function of the provided data rate level, where higher data rate levels would yield higher profit. Due to high complexity of the MINLP, we propose an efficient heuristic algorithm with lower computational complexity. We show that, when the UAV-BS location is determined, the resource allocation problem can be reduced to a Multidimensional Binary Knapsack Problem (MBKP), which can be solved in pseudo-polynomial time. To exploit this structure, the optimal bandwidth allocations are determined by solving several MBKPs in a search algorithm. We test the performance of our algorithm with two heuristics and with the MINLP model solved by a commercial solver. Our numerical results show that the proposed algorithm outperforms the alternative solution approaches and would be a promising tool to improve the total network profit.


I. INTRODUCTION
Unmanned Aerial Vehicle Base Stations (UAV-BSs) are expected to be used in the next generation Wireless Communications Networks (WCNs) for enhancing the capacity of the network as well as expanding the coverage [1]. Although the rapid deployment and mobility advantage of a UAV-BS has the potential to substantially improve the Quality-of-Service (QoS), there exist certain technical difficulties to be addressed The associate editor coordinating the review of this manuscript and approving it for publication was Fang Yang . such as resource management and channel modelling [2]. Among others, the location of the UAV-BSs play a key role to assist the WCN since positioning UAV-BSs optimally has vital importance when compared to terrestrial base station positioning. In traditional terrestrial networks, base stations are typically located with respect to long-term traffic estimations. Even though these estimations vary in time, installing new cells should be a cost-effective choice to meet the demand. On the other hand, a UAV-BS's location can be adjusted to meet instantaneous demand or to enhance the capacity of the WCN, yet, determining the optimal 3-D locations of such vehicle remains as a challenging and complex problem.
Not only does determining the locations of the UAV-BSs improve the network performance, but also resource management should be jointly addressed to boost the benefit of UAV-BSs. Similar to resource management strategies in classical wireless communications such as cache-enabled base stations [3], it is a significant problem in UAV-assisted WCNs to efficiently manage resources. Therefore, there has been a number of attempts to jointly optimize the location and resource allocation decisions [4]. Indeed, it is important to take the finite backhaul capacity of the UAV-BSs in consideration. Since the capacity is a significant concern for agility and reliability, joint optimization of the location and allocation decisions is imperative for achieving realistic assessment of the benefits of UAV-BSs.
In this study, we, specifically, address the aforementioned drawbacks and consider a UAV-assisted WCN with a capacitated UAV-BS serving to users on the ground who are not able to receive service from ground base stations (GBSs), e.g., due to high path loss. While serving the users, the backhaul capacity of the UAV-BS and the available bandwidth of the WCN are also considered in the problem setup. The objective is to maximize the total system profit. In particular, the users are assumed to be offered different data rate levels to improve their satisfaction and each user is allowed to select at most one of the offered options. This approach can also be used as a new pricing model in the next generation WCNs since the extent of the offered service type is a significant factor to ensure both customer loyalty and satisfaction [5].
Moreover, many studies related to the UAV-BS location have assumed that either the altitude or the projection of the utilized UAV-BSs on the ground are fixed. Such approaches are shown to transform the problem into a lower computational complexity case. However, such simplifications would cause sub-optimal decisions. For instance, it is shown that allowing the UAV-BSs to move both vertically and horizontally boosts the network performance in terms of throughput, resource utilization, and coverage [6]. Therefore, an agile approach should be developed in which the assumptions that leads to sub-optimality are avoided to have more realistic results.
Although there have been some attempts to consider network profit as an objective in 5G [4], this concept still needs to be thoroughly analyzed as the next generation WCNs are expected to embrace different pricing models. For instance, congestion-based [7] or tiered pricing models [8] are shown to be effective models to mitigate the burden on resource allocation decisions in the networks by charging higher prices for users whose demand are exhausting the network resources. These results motivate us to employ a service-level-based pricing model, where users are assumed to have different willingness values for different service levels. The service provider tends to serve users whose willingness values are higher, with a single service level. By adopting a servicelevel-based pricing model, the optimized resource allocations are expected to increase the profit and efficiency of resource management for service providers.
There exist different pricing models in the literature as can be seen in [9], [10] and the references therein. Some of these models include dynamic schemes, in which buyers, e.g., users, and sellers, e.g. service providers, sequentially make decisions on resource allocations and service provision, such as auction-based or priority-based pricing [10]. Some other models include static schemes, in which the prices of different services are set at the beginning of a finite planning horizon and the sellers provide the agreed service-levels to the buyers with no update on the predefined prices until the planning horizon terminates. In this study, a static service-level-based pricing scheme is employed since the proposed system model captures a single snapshot of a finite time horizon.
A straightforward application of the proposed model can be seen in scenarios where the exact values or at least the distribution of user willingness values is known a priori, and a service provider wishes to set the prices of different instantaneous data rate levels. In this case, the service provider would determine the optimal subset of users to serve under a limited backhaul capacity and set the prices of each data rate level with respect to this optimal subset. Since a user is typically assumed to accept a service if her/his willingness value is greater than the service price [11], the service provider can set the price of each data rate to the willingness value of the user whose willingness value for that specific data rate level is the minimum among the selected users, so that all users in the selected subset are guaranteed to be served.
Note that the aim of this article is not to compare how different pricing models perform. Instead, we provide an optimization problem for a single network operator who compete in a selfish market and offer a range of data rate levels to maximize its profit based on a static service-level-based pricing model [10]. Nevertheless, we include different static pricing models in the computational study and empirically show how such models can affect the network profit. We summarize the contributions of this article as follows: • Backhaul-aware optimization: Our proposed system considers both the backhaul and access capacities. Most of the studies related to the UAV-BS location have studied uncapacitated UAV-BS or assumed only one of these capacities, which would result in infeasible or unrealistic solutions. Therefore, we include both capacities in our problem setup so that the proposed model can be easily applied to real-life cases.
• Profit maximization: QoS-based optimization problems have covered a wide range of performance indicators such as throughput, latency or coverage [12]. Although these indicators seem reasonable, there has been little incentive for the network operators to increase the willingness to use UAV-BSs. We deliberately adopt the network profit as the objective function in the model so as to motivate network operators.
• Efficient solution approach: We use MINLP techniques to formulate the proposed model. Due to high complexity of this formulation, finding optimal solutions is difficult. We develop an efficient heuristic algorithm, which solves the problem within a reasonable time. We also provide the complexity of this algorithm.
• Computational study: We perform an extensive computational study in which we compare the performance of our algorithm with two more heuristics and one of the well-known MINLP solvers, the BARON solver over synthetically generated data. Our numerical results prove that our algorithm outperforms all three of the solution approaches. Moreover, we present empirical results to demonstrate the impact of pricing factors such as the cost of bandwidth, and user attitude in terms of willingness values, on various performance indicators from network profit to throughput.
• Theoretical contributions: We provide two important proofs regarding the concavity and unimodality of the data rate function, which is utilized extensively in the UAV-BS literature. More specifically, (i) we prove that the data rate function is concave with respect to bandwidth and (ii) we prove that the data rate function is unimodal with respect to the UAV-BS altitude. These results form the basis of proposed heuristic algorithm. The rest of the paper is organized as follows. Section II presents a literature overview on UAV-BS location problems. The system model and mathematical formulation of the proposed system are given in Section III. In Section IV, we propose the heuristic algorithm. Section V presents the computational results and Section VI concludes the paper with several future research directions. Proofs of concavity and unimodality of the data rate function are presented in Appendices A and B, respectively.

II. LITERATURE REVIEW
The opportunities brought by the UAV-BSs (such as enhancing capacity, improving QoS, extending coverage) are so promising that despite the relatively recent appearance of the topic, the literature on UAV-BSs has grown rapidly. Especially, the location optimization problems have attracted significant interest since they have a significant impact on the network performance [13]. In this section, we present a high-level overview of the UAV-BS location literature and present the differences of our study from the existing literature.
Recent studies on UAV-BS location have considered several design challenges such as 3-D deployment, air-to-ground channel modelling, and trajectory optimization. For instance, [14] jointly minimizes the number of UAV-BSs and determines the locations in areas with different user densities. In [15], a mathematical model is proposed for the single UAV-BS location problem, in which the ratio of the altitude of the UAV-BS to the radius of the area is used to reduce the problem to a 2-D location problem. In [16], the location of an uncapacitated UAV-BS is determined to satisfy different QoS requirements of the users. In [17], the location of multiple UAV-BSs is explored over a finite area to enhance the coverage and lifetime of the UAV-BSs. In [18], a joint optimization problem in which the locations of the multiple UAV-BSs and the association of the users to the UAV-BSs are determined while the network is assumed to store caching information of the users. In [19], the minimum average data rate provided to the users in a finite time horizon is maximized by jointly determining the trajectories of the UAV-BSs and the user-UAV-BS associations such that the maximum hovering distance is not violated. In [20], a deep reinforcement learning model is developed to maximize the system reward, which is found as the sum of sigmoid functions derived from difference between offloaded tasks and energy consumption of the UAV-BS. In fact, these studies, implicitly, assume that the backhaul capacity does not create a bottleneck when optimizing the UAV-BSs which is justifiable in some cases.
While there exists a number of studies that investigate backhaul-aware UAV-BS location optimization [14], [21], [22], the network profit aspect of the UAV-BSs have not been considered. In [21], the location of a single UAV-BS is considered to maximize the coverage such that the backhaul and access capacities are not exceeded. In [22], the locations of multiple UAV-BSs are determined to provide wireless service to the users who have a predefined delay-tolerance parameter. The resource allocation to the users in this study is also considered and an exhaustive search procedure is proposed to maximize a logarithmic utility function. However, in both studies, it is assumed that the users are identical in terms of QoS requirement. Since individual users are likely to have different requirement, a more realistic scenario would be assuming different demand values for users.
While the literature on UAV-BS location is rich, there has been surprisingly little work on investigating how pricing schemes can boost the QoS in these networks. It is shown in [23] through comprehensive user trials on the quality perception that users are likely to pay for enhanced network quality. Hence, paying more for better service can be considered as a typical customer tendency. Aligned with the model we propose in this study, it is shown that a well-designed pricing model, that does not depend only on flat-rate usage but also depends on the diversity and quality of services, can help network operators to improve profitability [24]. Pricing models can also help network operators improve the QoS. For instance, determining higher price levels for services that cause congestion could help network operators take control of excessive demand [25].
To date, various pricing models have been studied in wireless networks [9]. Dynamic and static pricing schemes have been used in various applications from user association [26], [27] to spectrum allocation [8], [28] and interference management [29]. We refer the interested readers to [9] and [10] and references therein for more detailed pricing models. We also refer the interested readers to Table 1 in [12] to see optimization problems that include pricing models with backhaul considerations in the problem setup.
Even though the existing studies have considered many different aspects of pricing strategies in WCNs, similar studies in the UAV-assisted WCNs, however, have been, mostly, overlooked in the literature. The related literature has been restricted to either uncapacitated network structures or fixed service-levels. In [4], a heterogeneous wireless network with UAV-BSs is designed to improve network profit. In [11], the UAV-BSs are considered to serve a number of given hotspots with random users and the optimal pricing for a single service-level is studied. However, both the backhaul and access capacities, and users' QoS requirements play significant roles for agile and reliable service provision. Therefore, we include both capacities and include varying QoS requirements in our system model.
The proposed model in this article covers the three pillars of a well-designed WCN, which has not been jointly considered to date: the UAV-BS deployment, backhaul and service-based pricing model. Since the state-of-the-art in UAV-assisted WCNs have mainly focused one or two of these pillars, the existing models and solution approaches cannot suit to a model in which all three pillars are combined. We, therefore, develop a novel UAV-assisted system model in the next sections and propose an efficient solution approach to the developed model.

III. THE SYSTEM MODEL
In our model, we consider a UAV-assisted WCN which consist of multiple GBSs and one UAV-BS. The UAV-BS has a dedicated backhaul link with one of the GBSs to serve the users on the ground. We assume that some of the users in the network cannot be served by GBSs due to several reasons (e.g., congestion in the wireless network or low channel quality). The UAV-BS is used as a temporary service provider for those users. The users are offered different data rate options and allowed to select at most one option. The profit is earned based on the selection of a specific user. Each user can have a different willingness value to pay for the same service. This approach can be considered as a variation of ''Paris Metro Pricing'' [9], in which the cumulative distribution function of willingness of the users is assumed to be known by the network operator instead of the exact willingness values. Fig. 1 illustrates a sample representation of the considered system with two data rate options: low and high. The willingness values of the users to pay for the service are shown by the number of $ icons in the figure. Although, the UAV-BS can rapidly be deployed anywhere in the air to cover as many users as possible, the availability of a backhaul link should be considered carefully to provide a reliable connection. We assume that each GBS has sufficient connection capacity with the communications infrastructure to transmit the data received from the UAV-BS to the core network. Therefore, there is no congestion in the fixed infrastructure.
The presented model can be validated based on the idea of enhancing the capacity as well as providing rapid supply to difficult-to-predict situations such as crowded events and activities. For instance, the area enclosed within the larger circle in Fig. 1 is a potential use case of such a model. The user distribution is not homogeneous, that is, some parts of the region hosts more users than the rest of the region. The UAV-BS can serve users in the congested areas more efficiently than GBSs. The presence of a UAV-BS also allows the network operators to utilize the idle capacity when the demand is not as high as the supply in some regions. Hence, the network operators can allocate the idle resources to the UAV-BS to create a relay connection between the lightly loaded GBS and the users who cannot be served directly by that GBS.
In the rest of the paper, I = {1, . . . , n}, J = {1, . . . , m} and K = {1, . . . , s} denote the set of users, GBSs, and offered data rate options, respectively. We assume that the service area of the UAV-BS, Q ⊆ R 3 , and the area in which users are located, S ⊆ Q, are finite. We use X ∈ Q to denote the location of the UAV-BS to be optimized. y u i ∈ S and y g j ∈ S denote the given locations of user i and GBS j, respectively. Moreover, the available bandwidth of GBS j that can be allocated to the UAV-BS is assumed to be b g j ≥ 0. The horizontal and vertical distances between two points, x, y ∈ Q, are found as r = ||L(x − y)|| and h = ||M (x − y)||, respectively, where L and M are linear transformations from R 3 to R 2 and R, respectively, and ||·|| is the l 2 norm. Throughout this article, unless otherwise specified, we use boldface capital letters to denote matrices and lower-case letters to denote vectors consisting of scalar parameters or variables denoted by the same letter (e.g., b g ∈ R m is the vector whose components are b g j ≥ 0 for j ∈ J and Y u ∈ R n×3 is the matrix whose components are y u i ∈ S for i ∈ I ). Moreover, the decision variables are denoted by upper-case letters, while the parameters are denoted with lower-case letters.

A. CHANNEL MODEL
There exist various air-to-ground channel models proposed in the literature [30]. However, the most widely employed model is presented in [31], therefore, we adopt this model. According to this model, there exist two propagation regimes such that the users in the first regime can have Line-of-Sight (LoS) connections and the users in the second regime have Non-Line-of-Sight (NLoS) connections and can maintain their connections due to the mechanisms of electromagnetic wave propagation which can be used to convey information beyond the obstructions (e.g., reflection and diffraction).
The probability of LoS is defined as a function of the locations of the UAV-BS and the users (i.e., P LoS : Q × S → [0, 1]). This probability for a specific user located at y u ∈ S when the UAV-BS is located at X ∈ Q can be calculated as where α and β are environment-specific constant parameters and θ (X , y u ) = (180/π ) arctan(h(X , y u )/r(X , y u )) is the elevation angle. Then, the pathloss between a specific user located at y u ∈ S and the UAV-BS located at X ∈ Q is defined as where d(X , y u ) = ||X − y u || is the distance between X ∈ Q and y u ∈ S, A = 10η log 10 (4π f c /c) + µ NLoS , and B = µ LoS − µ NLoS are constant parameters with f c denoting the carrier frequency in Hz, c denoting the speed of light in m/s, η denoting the path-loss component, µ LoS and µ NLoS denoting the associated excessive pathloss in dB with probabilities P LoS and P NLoS , respectively, and P NLoS = 1 − P LoS [31]. Note that the first term in the first equation in (2) denotes the free-space pathloss.
Unlike the users, the GBSs are generally located at carefully determined and advantageous locations. This enables them to communicate with the UAV-BSs through a much better channel. Therefore, the backhaul link is assumed to always have LoS communication with the GBSs, i.e., P LoS = 1. Then, the pathloss between the UAV-BS located at X ∈ Q and a GBS located at y g ∈ S is calculated as L g (X , y g ) = A + 10η log 10 d(X , y g ) + B. (3)

B. PROBLEM FORMULATION
We consider point-to-point wireless connections between the UAV-BS and the GBSs. This connection is assumed not to interfere with user links of the UAV-BS. To facilitate such an assumption, reversed time-division duplexing is employed to avoid interference between the backhaul and user links, such that, during downlink of GBS, UAV-BS is in the uplink mode. However, when GBS is in the uplink mode, users may largely be affected by the GBS due to the self-interference. To avoid this, orthogonal frequency channels are used in the backhaul and user links [32]. The objective in the proposed model is to maximize the network profit by jointly determining the UAV-BS location and allocation of the available bandwidth to the users subject to the backhaul and access capacities. We first give the definitions of the data rate and backhaul capacity and then introduce the associated Mixed-Integer Non-Linear Programming (MINLP) formulation.
For a specific user located at y u ∈ S, the actual data rate, when the allocated bandwidth is equal to B u ∈ [0, B max ] (in Hz) from the UAV-BS located at X ∈ Q, where B max is the maximum bandwidth amount that can be allocated to a user in practice, can be found as where with p d and ω N denoting the transmit power of the UAV-BS and a noise figure, respectively. Note that the data rate function, R, is concave with respect to B u (see Appendix A for the proof) and unimodal with respect to the altitude of the UAV-BS (see Appendix B for the proof). This geometry will become a significant component of the solution approach we propose in Section IV.
The backhaul capacity of a UAV-BS located at X ∈ Q, when it has a backhaul link with a GBS located at y g ∈ S, can be found as where S g (X , y g ) = [p g − L g (X , y g ) − 10 log 10 b g − ω N ]/10 is the SNR of the UAV-BS with p g denoting the transmit power of the GBS. Recall that the available bandwidth at GBSs, b g , is assumed to be known a priori.
Slightly abusing the notation, we use g j ) to denote the actual data rate of user i ∈ I and the backhaul capacity provided to UAV-BS from GBS j ∈ J when the UAV-BS is located at X ∈ Q, respectively. We introduce the binary variable, Z j ∈ {0, 1}, to indicate whether or not the UAV-BS has a backhaul link with GBS j, and the piece-wise linear function, u i (R i ) : R → R, to denote the profit gained when user i ∈ I is served with data rate R i ≥ 0. Consequently, the MINLP formulation is defined as The objective function in P maximizes the total network profit, constraints (6) and (7) ensure that the backhaul and access capacities are not exceeded, respectively, while constraint (8) requires the UAV-BS to have a backhaul link with exactly one GBS.

C. PRICING MODEL
As explained in Section I, we adopt a Paris Metro Pricing model, in which users are assumed to have different willingness values to pay for different data rate options. Therefore, VOLUME 8, 2020 the profit gained from user i ∈ I is modeled as where δ k denotes the k th data rate option in K , and φ ik denotes the willingness value of user i ∈ I to pay for the data rate option k ∈ K . The φ values are assumed to be an ascending order (i.e., φ i1 ≤ φ i2 ≤ . . . ≤ φ is ), thus, providing higher data rate to a specific user would yield higher profit. This approach is based on the fact that the price of a service, generally, increases when its quality is improved [23]. Moreover, each user is assumed to have different willingness values for the same option, so that a diverse population can be represented.
Note that each user has a single bandwidth amount to be served with a specific data rate since the data rate function is concave in B u . As K has a finite number of data rate options, i.e., s < ∞, and the u function has a piece-wise structure, where the profit gained from a user does not change for a finite range of data rate values, allocating more bandwidth to a user after satisfying a specific data rate option in K does not yield additional profit until the next data rate option is satisfied. For instance, there would be no difference in the profit between allocating more bandwidth to provide 1.1 Mbps instead of 1 Mbps to a user unless the u function has different values at 1 Mbps and 1.1 Mbps data rates. Since the objective of (P) is to maximize the total network profit, bandwidth allocations are determined by providing the least amount of bandwidth to users for each data rate option so that remaining bandwidth can be allocated to some other users to improve network profit. Therefore, the bandwidth allocations are upper bounded, implicitly, in (P).
P is a complex problem to solve since it includes a non-convex constraint set as well as binary and continuous variables. In fact, it belongs to the NP-complete problem class since relaxing the backhaul capacity constraint (6) and fixing the bandwidth allocations would reduce the problem to the maximal covering location problem which is known to be NP-hard [33]. Therefore, in the next section, we propose an iterative solution approach with lower complexity.

IV. SOLUTION APPROACH
In this section, we propose a heuristic algorithm in which P is solved in an iterative manner. Note that the backhaul capacity can be explicitly found when the UAV-BS location is fixed, which substantially decreases the complexity of the problem. We exploit this relaxation and develop an efficient algorithm in which the altitude and the horizontal coordinates of the UAV-BS are searched and at each fixed altitude several resource allocation problems are optimally solved for some fixed coordinates.

A. PRELIMINARIES
Let C j (X ) be the actual backhaul capacity received from GBS j given the UAV-BS location, X ∈ Q. We introduce the binary variable, T ik ∈ {0, 1}, to indicate whether or not user i ∈ I is served with data rate option k ∈ K . There does not exist a closed form expression to determine the required bandwidth to serve a user located at y u ∈ S for a fixed data rate value, δ. However, since we proved that R is concave with respect to B u , this bandwidth can be explicitly found with a line search algorithm (e.g., bisection search). Let B u ik (X ) denote the required bandwidth to serve user i ∈ I with data rate option k ∈ K when the UAV-BS location is given as X ∈ Q. Then, P can be reduced to the following Integer Linear Programming (ILP) formulation, The objective function in P maximizes the total network profit for a fixed UAV-BS location. Constraints (10) and (11) ensure that the backhaul and access capacities are not exceeded. Constraints (12) stipulate that each user can only be served with at most one data rate option, while constraint (13) guarantees that the UAV-BS has a backhaul link with exactly one GBS.
The most important advantage of P is that the formulation is now free of the non-convex data rate and backhaul capacity functions. In this way, instead of determining the actual data rate values of users with the location and allocation variables, the problem is altered to an assignment problem in which the data rate options are assigned to users regarding their profit values. Note that, P can be solved separately for each GBS j ∈ J . The optimal solution is found when the UAV-BS has a backhaul link with the GBS that provides the highest backhaul capacity, i.e., the highest value in the right-hand side in (10). This would allow the UAV-BS to serve more users in case all other variables are fixed. Therefore, it is sufficient to solve this relaxed problem by setting the Z variable corresponding to the GBS that yields the highest backhaul capacity for the given UAV-BS location to 1 and all other Z variables to 0.
Let j * (X ) be the GBS providing the highest backhaul capacity for a given X , i.e., j * (X ) = arg max j∈J C(X , y g j , b g j ). Then, omitting the Z variables and associated constraint (13), P can be reformulated as i∈I k∈K The resulting P j * (X ) is also an ILP formulation, which can be efficiently solved by commercial solvers such as CPLEX or GUROBI. In fact, this problem is known as ''Multidimensional Binary Knapsack Problem'' (MBKP), which can be solved in pseudo-polynomial time via dynamic programming techniques [34]. In MBKP, the problem is to select a subset of items subject to multiple capacity constraints. The objective is to maximize the total utility of the selected items that have specific utility values such that none of the capacity constraints are exceeded. Regarding our formulation in P j * (X ) , items can be considered as user-data rate pairs, i.e., T variables, while there exist two types of capacities, i.e., constraints (14) and (15). Therefore, items can be selected until either backhaul or access capacity is fulfilled.
By exploiting the structure explained above and our proof of the unimodality of the data rate function with respect to the UAV-BS altitude, we propose an iterative heuristic algorithm in the following subsection.

B. HEURISTIC ALGORITHM
We proved that the data rate function is unimodal with respect to the UAV-BS altitude. As a consequence of this fact, our heuristic algorithm uses a line search over the range of altitudes. At each altitude level, another search algorithm is used in which the revised relaxed formulation, P, is used to determine the optimal bandwidth allocations. The summary of the algorithm is given in Algorithm 1.
We utilize the well-known ''Golden Section Search'' (GSS) algorithm to search the optimal altitude of the UAV-BS. Before explaining the details of the algorithm, we provide a brief description of the GSS. The GSS algorithm is typically used to find the global minimum or maximum of a unimodal function by iteratively searching the function domain. Suppose we are given a function, f : R → R, that is assumed to have a single maximum within [q, w], where q, w ∈ R and −∞ < q < w < +∞. Then, two interior points are selected by using the golden ratio, γ = √ 5−1 2 . The first point, x 1 , is equal to q + γ (q − w), and the second point, x 2 , is equal to w − γ (q − w). If f (x 1 ) ≥ f (x 2 ), then it is sufficient to say that the global maximum is between x 2 and w and the same procedure continues within [x 2 , w], otherwise, the search continues within [q, x 1 ]. This recursive procedure continues until the range size drops below a predefined approximation threshold, g . As a result, the center of the final range is determined to be the maximum of f . end if 10: end while 11:

Algorithm 1 Finds the UAV-BS Location and the Bandwidth Allocations of the Users
. 13: return X * , B * .
Based on the above procedure, our algorithm determines two solutions at each iteration for the original problem by fixing the UAV-BS altitude to two different levels. Let h l and h u denote the minimum and maximum altitudes to which the UAV-BS can be located, respectively. The two altitude levels are determined as h 1 and h 2 , i.e., For each altitude level, a grid search algorithm, which is summarized in Algorithm 2, is applied to determine the best UAV-BS location and optimal bandwidth allocations at the corresponding altitude (Lines 3-4 in Algorithm 1). Let B 1 and B 2 , respectively, denote the optimal allocation decisions at X 1 and X 2 , which are the best locations at altitudes h 1 and h 2 , respectively. These locations and decisions are used to determine the direction of the GSS. If X 1 yields a higher profit than X 2 with the corresponding bandwidth allocations, then the altitude range is updated as [h 2 , h u ], otherwise, this range is updated to [h l , h 1 ]. The GSS algorithm terminates whenever the size of this range drops below g . Note that as a last step, the relaxed problem is solved once more, where the UAV-BS altitude is equal to the center of the final altitude range (Lines 11-12 in Algorithm 1), and the best solution found at this altitude is reported.
In Algorithm 2, for a fixed altitude level, h, a grid set, S(h), is generated by dividing S into Q smaller grids, Algorithm 2 Searches the Horizontal Plane for a Fixed UAV-BS Altitude Input: Y u , Y g , K , φ, S(h).
1: fors q ∈ S(h) do 2: j * = arg max j∈J C(s q , y g j , b g j ). 3: Solve P j * (s q ). Let T * q be the optimal solution. Set q ← (T * q |s q ). 4: end for 5: q * ← arg max q=1,...,Q q . 6: returns q * , B u (q * ). VOLUME 8, 2020 i.e., S(h) = {s q ⊆ S : h = h, q = 1, . . . , Q}, with s q denoting the center of grid q. For each grid, the relaxed problem, P, is solved, where the UAV-BS location is assumed as the center of the corresponding grid, i.e., X =s q , q = 1, . . . , Q. Then, the objective function values attained from all grids are compared and the center of the grid that yields the highest profit is returned as the best objective function value for h. In this algorithm, q denote the objective function value attained from grid s q .

C. CONVERGENCE AND COMPUTATIONAL COMPLEXITY
The developed algorithm terminates after a finite number of iterations since the altitude range is assumed to be finite. Suppose that the initial altitude range is [h init l , h init u ]. Then, after each iteration, this range gets narrower by a constant rate, γ ≈ 0.618, thus, after k iterations, the altitude range is decreased to (h init u − h init l ) × γ k . As k goes to infinity, this range goes to 0. Since g > 0, the algorithm converges after a finite number of iterations.
From the computational perspective, we solve several MBKPs in each iteration, whose complexity depends on the number of items and the capacities [34]. As we have n × s items for each problem and two capacities, C and b g , the complexity of solving a single MBKP is O(ns(C + 1) (b g + 1)). The number of MBKPs to be solved for each fixed altitude is determined by Q. Note that providing a dense grid structure would improve the solution accuracy, while it may substantially decrease the computational performance and vice versa. As a result, the overall complexity of our algorithm is O(R − Qns(C + 1)(b g + 1)), where R − is the number of iterations to achieve g approximation, i.e., R − = min{R ∈ Z + : (h init u − h init l )γ R ≤ g }, where Z + denotes the non-negative integers.

V. COMPUTATIONAL RESULTS
In this section, a detailed computational study is performed to evaluate the performance of the developed algorithm in terms of solution accuracy and computational efficiency. We also provide some powerful insights on how service-based pricing can boost the network profit. All simulations are performed for a suburban area of 1500 m×1500 m on an Intel i7-6700 CPU @3.40 GHz, 64-bit, 8GB RAM Windows 10 computer, where the communication and other simulation parameters are given in Table 1 as suggested in [14] and [22] for urban environment.

A. BENCHMARKING TOOLS
For comparison, we use the optimal solutions (or the best feasible) attained from the well-known MINLP solver, the BARON solver [35], provided on the NEOS platform [36] with 1-hour time limit and two other heuristic approaches. Both our algorithm and heuristic algorithms are coded in Python v3.6.
For the first heuristic algorithm, namely ''Heuristic R '', the UAV-BS location is randomly determined in Q and optimal bandwidth allocations are determined by solving P at this random location. To smooth the randomization effect, we use the best objective function values attained after 50 replications. For the second heuristic algorithm, namely ''Heuristic F '', instead of a random location, the UAV-BS is located to a fixed point, where the horizontal coordinates are determined as the weighted average coordinates of user locations and the altitude is determined as the quarters of the altitude range yielding four different altitudes. Consequently, the optimal bandwidth allocations are determined by solving P at these four fixed locations and the location yielding the highest profit value is selected as the UAV-BS location.

B. DATA GENERATITON
We assume that the GBSs are uniformly located in S and the UAV-BS altitude range is between 50 and 500 meters. The users are assumed to be located according to a Poisson Point Process (PPP). To apply PPP, the users are randomly located to specific parts of the region on S with a random number of parent nodes and clustering rate. For instance, if the number of parent nodes is 5 and the clustering rate is 70% for an instance with 100 users, then 70 users are located in close proximity of uniformly selected 5 different points while the remaining 30 users are located uniformly in S. The offered data rate options are given in Table 2. The willingness values of the users are determined in a way that higher data rate values yield higher profit, i.e.,   all the random numbers in the replications are generated by scipy.random libraries. Fig. 2 depicts the change in the normalized objective function values for all solvers with respect to different number of users for each data rate option set. The average CPU time of our algorithm is 255.1 seconds, while it decreases to 29.4 and 6.3 seconds for Heuristic R and Heuristic F , respectively. The BARON solver reports the optimal solutions for only instances with 50 users and 2 data rate options, thus, we can only report the objective function values of the best feasible solutions whenever possible. Fig. 2 clearly shows that our algorithm outperforms all other solvers in terms of total profit. The Heuristic R is the worst performing algorithm because of the insufficient tolerance of avoiding local optima in determining the UAV-BS location. In addition, the BARON solver suffers from increasing number of binary variables. Since the BARON solver uses a branch-and-bound technique over the binary variables, increasing number of these variables would significantly increase the number of branches and it would be inefficient to solve resulting non-convex optimization problem even after fixing some of the binary variable values on each branch. Note that our algorithm finds the optimal solutions for those instances for which the BARON also finds the optimum (2 instances and 20 replications).

C. PERFORMANCE ANALYSIS
An interesting result observed from Fig. 2 is that the Heuristic F performance is better than that of Heuristic R although it is outperformed by our algorithm and BARON. The main drawback of Heuristic F is that it is not always the best alternative to locate the UAV-BS to a central point, instead, optimally determining the location aligned with bandwidth allocations leads to a higher network profit. Nevertheless, the CPU times of Heuristic F is lower than the others which suggests that it can be employed to find a fast primal solution for large-scale problems.
Heuristic F performance worsens with the number of data rate options. This is due to the trade-off between allocating more bandwidth to the users who yield higher profit in return and following a fairer allocation where the available bandwidth is allocated according to the central position of the UAV-BS. Since the objective in our case is the network profit, the latter is advantageous. We have also observed that Heuristic F performance further worsens with higher clustering rates and parent node numbers since the weighted average location is unlikely to represent the overall dispersion of the users in such cases. Therefore, exploring the horizontal plane (e.g., grid search) appears to be a reasonable approach to improve the network profit.
In Fig. 3, we present the change in coverage and total profit attained with respect to different data rate options. Note that the coverage is defined as the ratio of the total data rate provided to covered users over the maximum total data rate (i.e., n × δ s ). Fig. 3 shows that there exists a significant trade-off between improving the coverage and increasing network profit. As the likelihood of having more users with higher willingness values increases with more users, our algorithm prefers allocating the available bandwidth to those users who bring more profit which, in turn, leads to lower coverage.   Fig. 4 illustrates the change in the total profit between the solution found by our algorithm for the original problem in which more than one option is offered and a variant of the original problem where users are assumed to be offered a single data rate (that is equal to the average data rate offered in the original problem) and the profit gained from a specific user is assumed to be equal to the average willingness value of the corresponding user. The figure implies that the improvement level would increase with increasing number of users. For instance, offering three data rate options to 50 users (e.g., 1, 2 and 4 Mbps) instead of offering a single option (e.g., 2.33 Mbps) improves the total profit by 44.6%. This value escalates to 49.9% when the number of users increases to 200. Note that this outcome is valid for the assumption proposed in our system model that users are willing to pay more for improved service. Nevertheless, this assumption is highly likely to occur in real-life cases since a customer does not consider only the price but also does pay attention to the service quality [5], [23].

D. PRICING ANALYSIS
In this section, instead of comparing the performance of different solution approaches, we present insights on how different willingness values and pricing models affect the network performance. First, we vary the willingness values of users by incorporating different probabilistic factors into the procedure that we have already applied in the preceding analyses. Second, we revise the pricing model by adding a new cost factor corresponding to the allocated bandwidth. In particular, we modify the utility function defined in Eq. (9) such that the allocated bandwidth has a negative impact on the profit gained from a specific user. We then resolve the same instances that are used in the preceding analyses and present the results. Since our algorithm is shown to outperform all other approaches, the comparison is based on only results obtained by our algorithm.

1) EFFECT OF WILLINGNESS VALUES
Note that all the previous analyses are based on the assumption that willingness values of users are determined with a piece-wise increasing function since it is assumed that higher willingness values can be perceived as a tendency to the better service-level. We keep this assumption but with a difference, where, instead of using only uniform distribution, we consider the normal, log-normal, and gamma distributions to generate random numbers in the willingness values of users.
It is suggested in [9] that Paris Metro Pricing can be simulated with known probability distributions of willingness values although the exact values are not known. Hence, we employ the normal distribution to represent a user population with no aggressive differences, while the log-normal distribution allows us to represent aggressive marginal values. The gamma distribution is, on the other hand, selected to allow us to combine both approaches. Table 3 shows how different distributions are included in data generation, where N (µ, σ ) and L(µ, σ ) are normal and log-normal distributions with mean µ and standard deviation σ , respectively, while G(α, β) is the gamma distribution with shape parameter α and rate parameter β. All other parameters such as user locations and pathloss parameters remain unaltered. We also include a base scenario in our analysis, where all users are assumed to have a single willingness value regardless of the data rates offered. For this case, each instance generated in Section V-C is modified as if each user is associated with the average of willingness values. Fig. 5 illustrates the change in the coverage, profit, and data rate per allocated bandwidth amounts with respect to different distributions for each number of users. Since, especially, the random variables generated by the log-normal distribution are significantly greater than those generated by the normal and gamma distributions, all figures are presented after normalization. Fig. 5 clearly shows that log-normal distribution yields the highest profit values, while the base scenario results in the highest average data rate values. Coverage values follow the same trend for all distributions with slight differences which do not appear to be statistically significant. This result shows that unifying willingness values result in increasing coverage and average data rate values which confirms that there is no additional revenue in allocating more bandwidth to a specific user when no marginal profit is expected, hence, the access capacity is allocated more uniformly. Moreover, offering different service levels to a population which is unlikely to pay more for additional data rate does not seem to be a fruitful approach. Instead of offering new services in such case, investing in the access capacity appears to be promising in serving more users.
An interesting result from our analysis is that when users are expected to have diverse willingness values as it is represented in the log-normal case, offering new services with higher quality is preferable. Although increasing trend in the network profit is obvious with all distributions (see Fig. 4 and Fig. 5b), the log-normal distribution results in the highest improvement since the marginal difference is expected to be higher.

2) EFFECTS OF INCORPORATING THE COST OF BANDWIDTH TO THE PROFIT MODEL
It is assumed in Eq. (9) that the utility gained from a user is defined as a function of data rate, only. In the following, we relax this assumption and add a cost factor to profit model. There exist a number of methods in the literature, in which the profit (or utility in our case) function includes the costs related to resource utilization, i.e., the resource allocated to users, such as bandwidth [4] or energy [37]. Since our model requires the available bandwidth to be allocated, we modify the original utility function by adding the cost of bandwidth. In particular, u function is modified as where λ is the unit cost of bandwidth. Observe that Eq. (16) is a generalization of Eq. (9) since Eq. (16) can be reduced to Eq. (9) when λ = 0. Note that the objective function of both P and corresponding sub-problems should also be accordingly modified. However, our algorithm does not need a major update since fixing the UAV-BS location in each iteration, i.e., X , would allow us to determine explicitly the required bandwidth amounts, B u ik (X ), for each data rate option k ∈ K to serve user i ∈ I , and include these values as parameters in the sub-problems. In particular, it is sufficient to modify only the objective function in P j * (X ) as We resolve the instances generated in Section V-C with uniform willingness values for different values of λ, i.e., λ ∈ {10 −6 , 20 −6 , 30 −6 }, to analyze the effects of bandwidth cost on the network performance. Fig. 6 shows the change in the network profit, throughput, i.e., total data rate provided, number of served users and bandwidth utilization for different λ values with respect to number of users. Note that λ = 0 in the figures depicts the results that have been found previously.
It is clear that λ has a shifting effect (i.e., shifting the profit down) on the network profit as illustrated in Fig. 6a which shows that the network profit decreases with increasing unit cost values regardless of the number of users and the data rates offered. However, this shift is not observed in terms of throughput. Fig. 6b shows that offering more data rate options decreases the negative impact of bandwidth cost in terms of throughput as it can be observed that the variation of throughput values reduces with increasing data rate options for each number of users value among different λ values.
A significant impact of the new profit model can be observed in terms of resource utilization. Fig. 6c and 6d together show that higher bandwidth costs have a negative impact on the total number of served users as well as on the bandwidth utilization. Note that the utilization values are the ratio of total allocated bandwidth to the available bandwidth. In fact, the algorithm prefers to allocate bandwidth to users who are likely to yield higher profit and prevents allocating excess capacity to those whose willingness do not satisfy the incurred cost. Offering higher data rate options is preferable in such cases as it encourages the users who can upgrade their service-levels, which in return create opportunities for other users who would not be served otherwise. Such an outcome can be observed in Fig. 6d that average utilization increases from 85.3% to 92.8% when 4 and 8 Mbps options are included in the data rate options.

VI. CONCLUSION
UAV-BS location problems have attracted significant interest in both industry and academia due to their potential to bring unprecedented advantages like rapid deployment, dynamic coverage extension, and on-demand capacity increase. Most of the studies to date have focused on improving QoS while the backhaul aspect and new pricing models are, mostly, overlooked. In this study, we address both the capacity aspect of QoS and a new pricing model based on different service levels. A novel mathematical model is developed and an efficient two-phase solution procedure, that combines GSS and grid search, is proposed for jointly optimizing the UAV-BS location and bandwidth allocations to the users to maximize the network profit. The proposed algorithm is shown to be capable of significantly improving the network profit.
This study can be extended by considering the temporal dimension of the problem since users typically relocate in time and the willingness values of the users, possibly, change at different time intervals. Such dynamics can lead to changes in the location, the hovering time, and the trajectory of the UAV-BS which necessitates the design of dynamic pricing policies. Another extension can be to combine energy-related performance metrics with the backhaul capacity since energyefficient operation of UAV-BS is required to increase service time of the UAV-assisted WCNs. Covering all users with the minimum number of UAV-BSs by considering new pricing models and backhaul capacity is also an interesting future research avenue.

APPENDIX A PROOF OF CONCAVITY OF DATA RATE FUNCTION WITH RESPECT TO BANDWIDTH
In this appendix, we prove that the data rate function is concave with respect to b i . By definition, a function is concave if the second derivative of the function is non-negative.
Recall that the data rate function is defined as where S u (X , y u , B u ) = [p d − L u (X , y u ) − 10 log 10 B u − ω N ]/ 10. Since we will prove the concavity for only bandwidth, B u , we can assume that the other variables are fixed, i.e., X = X , y u = y. After a number of algebraic manipulations, we have where = 10 p d −L u (X ,y)−ω N −1 > 0 is a fixed term. Having taken the second derivative of (17) with respect to B u , we have Since B u > 0, R (B u ) < 0 always holds.

APPENDIX B PROOF OF UNIMODALITY OF DATA RATE FUNCTION
In this appendix, we prove that the data rate function is unimodal with respect to the UAV-BS altitude. By definition, a function, f : R → R, is unimodal if there exists an m ∈ R such that for some value, δ ∈ R, it is monotonically increasing (decreasing) for δ ≤ m and monotonically decreasing (increasing) for δ ≥ m. In this case, the maximum (minimum) value of f (·) is f (m) and there are no other local maxima (minima). If the function is differentiable, proving that the derivative of the function is equal to 0 at a single point in its domain proves its unimodality. Here we will make use of this property.
Recall the data rate function is defined as R(X , y u , B u ) = B u log 2 1 + 10 S u (X ,y u ,B u ) , where S u (X , y u , B u ) = [p d − L u (X , y u ) − 10 log 10 B u − ω N ]/ 10. Since we will prove the unimodality for only UAV-BS altitude, h, we can assume that the other variables are fixed, i.e., y u = y, B u = B, r(X , y u ) = r. After a number of algebraic manipulations, we have where where we denote a partial derivative with '' '', e.g., f (h) = ∂f (h) ∂h . Note that the denominator in (20) is nonnegative, since f (·) is a product of two non-negative values plus 1. For a given bandwidth of B ≥ 0, this derivative can only be 0 if f (h) = 0. The first derivative of f (h) can be found as Since we assume h > 0, the first term in (21) is also greater than 0. Therefore, f (h) = 0 holds if the expression in brackets is 0. Before proceeding, we first derive m (h) as m (h) = K log(N )N g(h) g (h). (22) By substituting (22) into the expression in the brackets in (21), we have The last transformation in (23) follows since both K and N are greater than 0. Before proceeding, we first derive g (h) as where θ (h) = 180r π r 2 +h 2 . By substituting θ (h) into (24) and then (24) to (23) where A = 180(µ NLoS −µ LoS )rαβp ηπ ≥ 0, = αp ≥ 0, and q(h) = s θ(h) ∈ (0, 1]. For simplicity we use Q(h) = q(h)(A − 2 h) and W (h) = h 1 + ( q(h)) 2 to denote the left and right-hand side in (25). Since θ (h) is a monotonically increasing function in h and β > 0, q(h) = s θ(h) = e −βθ(h) is a monotonically decreasing function in h. When h starts to increase from 0, Q(h) monotonically decreases from A to 0. On the other hand, W (h) is a monotonically increasing function which increases from 0 to the maximum altitude, h + , when h increases. This implies that there can be found at most one h = h where Q h = W h which proves the unimodality of R with respect to h.
Also it is straightforward to show that R(·) monotonically increases for h < h and monotonically decreases for h > h, which implies that R(·) has its maximum value when h = h.
BULENT TAVLI (Senior Member, IEEE) received the B.Sc. degree in electrical and electronics engineering from Middle East Technical University, Ankara, Turkey, in 1996, and the M.Sc. and Ph.D. degrees in electrical and computer engineering from the University of Rochester, Rochester, NY, USA, in 2002 and 2005, respectively. He is currently a Professor with the Department of Electrical and Electronics Engineering, TOBB University of Economics and Technology, Ankara. His research interests include wireless communications, networks, optimization, machine learning, embedded systems, information security and privacy, smart grids, and blockchain.
HALIM YANIKOMEROGLU (Fellow, IEEE) received the B.Sc. degree in electrical and electronics engineering from Middle East Technical University, Ankara, Turkey, in 1990, and the M.A.Sc. degree in electrical engineering and the Ph.D. degree in electrical and computer engineering from the University of Toronto, Canada, in 1992 and 1998, respectively. He is currently a Professor with the Department of Systems and Computer Engineering, Carleton University, Ottawa, Canada. His collaborative research with industry has resulted in 37 granted patents. His research interest includes 5G/5G+ wireless networks. He is a Fellow of the Engineering Institute of Canada (EIC) and the Canadian Academy of Engineering (CAE). He is also a Distinguished Speaker of the IEEE Communications Society and the IEEE Vehicular Technology Society.