Interference coordination and throughput maximisation in an unmanned aerial vehicle-assisted cellular: User association and three-dimensional trajectory optimisation

Unmanned aerial vehicles in wireless communications have caught great attention in unmanned aerial vehicle assisted communications and cellular-connected unmanned aerial vehicles. However, the unmanned aerial vehicle’s advantage of strong line-of-sight channels also brings about severe and extensive interference between unmanned aerial vehicle and the terrestrial communication networks and between multiple unmanned aerial vehicles. This paper jointly optimises the unmanned aerial vehicle’s three-dimensional trajectory and ground users’ association schemes with the unmanned aerial vehicle and ground base station to reduce interference. In consideration of fairness, the authors formulate the optimisation problem as maximising the minimum average throughput among all ground users, a mixed-integer non-convex problem that is challenging to tackle. To ease the difﬁculties, the authors have the optimisation problem relaxed and propose a suboptimal solution in the framework of the block successive upper-bound minimisation algorithm. Further-more, an improved block successive upper-bound minimisation based penalty algorithm is obtained to avoid the local maximum point problem. Numerical simulations have validated our schemes.


INTRODUCTION
Unmanned aerial vehicles (UAVs) have recently been found promising for wireless communications due to their advantages of strong line-of-sight (LOS) channels, high mobility, and flexible deployment [1,2]. On the one hand, many types of UAVassisted communications have been proposed, such as aerial base stations (BSs) [3], UAV-enabled relays [4,5], multi-casting [6], and data collection [7]. On the other hand, UAVs can be used as aerial nodes to execute particular tasks while controlled by the terrestrial network called cellular-connected UAVs [8].
For example, a new kind of communication network named UAV-to-everything communications is put forward in depth [9, Moreover, one UAV could achieve multiple complex functions, as in [13], the UAV jointly optimises and performs sensing tasks and data transmission to the BS. However, the strong LOS links also cause severe and extensive interference, leading to widespread performance degradation and design complication [14,15]. Fortunately, a variety of researches have been carried out for interference coordination. For cellular-connected UAVs, the UAVs uplink cell associations and power allocations are jointly optimised in [14] to obtain optimal inter-cell interference coordination design and airground performance trade-off. A novel cooperative interference cancellation strategy by exploiting the existing backhaul connections is proposed in [16]. Moreover, a deep reinforcement learning algorithm is raised in [17] to design an interferenceaware path planning scheme for a network of cellular-connected UAVs, which is further discussed in a cellular Internet of UAVs in [18,19]. In [20], the author considers a cooperative cellular Internet that comprises multiple UAVs and has co-channel interference reduced through radio resource management and three-dimensional (3D) beamforming techniques. As for UAVassisted communications, a dynamic mobility-aware interference avoidance scheme is proposed by utilising spectral graph theory [21]. The ground nodes' transmit power allocation and scheduling and the UAV's 3D trajectory design are jointly optimised with the existence of a group of ground jammers in [22]. A UAV-enabled interference channel is considered in [23], where a joint trajectory and power control problem is formulated to maximise the aggregate sum rate. In [24], directional antennas at the UAV and ground base station (GBS) are jointly applied to eliminate the mutual interference.
To mitigate co-channel interference with low costs, we propose a method of joint optimisation without consuming extra communication resources. The GBS and UAV user association schemes are jointly designed with the UAV's trajectory, which enjoys a large degree of freedom. In this way, the UAV could always establish links with those ground users (GUs) less affected by the GBS. Furthermore, it would keep a distance from the GU that is communicating with the GBS at the same time. In respect of the GBS, it would switch to those GUs far away from the UAV in real time. With an appropriately designed algorithm, the targets of avoiding interference, improving throughput, and maintaining the fairness between users could be balanced and achieved in the meantime.
This paper investigates a UAV-assisted hybrid cellular network, where a UAV and a GBS are concurrently deployed to serve several GUs. The UAV and GBS utilise the same frequency band, so co-channel interference exists when they are simultaneously working. A time-division multiple access (TDMA) scheme is adopted. Having considered the fairness between GUs, a mixed-integer non-convex optimisation problem is formulated to maximise the minimum average throughput among all GUs under certain UAV trajectory and user association constraints. The UAV and GBS user association schemes, along with UAV's trajectory, are jointly optimised to coordinate the interference and maximise the throughput. Intuitively, the GBS would choose those GUs farther from UAV to reduce UAV's interference on its users, while the UAV would behave similarly. The formulated problem is quite challenging to solve, so we first have the integral constraints relaxed to be continuous. Then a suboptimal solution to the relaxed problem is proposed applying an algorithm in the framework of the block successive upper-bound minimisation (BSUM) algorithm [25]. Furthermore, an improved BSUM-based penalty algorithm is put forward to avoid the local maximum points problem. Our scheme and algorithms are validated to be effective by numerical simulations.
The rest of this paper is organised as follows. The background and principal work of our research are introduced in Section 1. In Section 2, we present the system model and formulate the optimisation problem to be solved. Two algorithms to tackle the formulated problem are proposed in Section 3, while the former is based on the BSUM algorithm, and the latter Information Signal Interference Information Signal Interference is improved with penalty factors. The convergence behaviours of these algorithms are also discussed. Numerical simulation results are provided in Section 4. Finally, this paper is concluded in Section 5.

System model
As shown in Figure 1, we consider a UAV-assisted single-cell wireless communication network. In the system, a UAV and GBS are jointly deployed to serve a group of K GUs. The GUs are generated by the mixed traffic model in [26]. Two sets of GUs are generated via an independent spatial Poisson point process with different densities u and nu . Then GUs with density nu are shifted towards GBS by correlation parameter p, meaning the distances between these GUs and the GBS are resized from R to R(1 − p). All GUs are indexed by the set  = {1, … , K }. A 3D Cartesian coordinate system is defined, where the GBS is fixed at q G = (x G , y G , H G ), and GU k is denoted by q k = (x k , y k , 0). The UAV is assumed to fly periodically and its position at time instant t is q(t ) = (x(t ), y(t ), z(t )). For ease of exposition, all the x and y coordinates are combined to form 2D vectors w = (x, y).
To avoid ground obstacles and obey government regulations, the UAV's altitude is restricted. The UAV's horizontal and vertical velocities are assumed to be controlled independently, and they cannot exceed their maximum limits V h max and V v max . Thus, we obtain the following constraints on UAV's trajectory: where T represents the flight period, and (1) is the periodicity requirement. T is equally discretised into N t time slots for ease of design, where the element slot length t is chosen to be sufficiently small that the communication conditions can be considered approximately unchanged during it. Therefore, the UAV's trajectory q(t ) can be approximated by a sequence {q[n]}, where q[n] = q((n − 1) t ) represents UAV's position at time slot n (n ∈  = {1, … , N t + 1}). Then the trajectory constraints (1)-(4) can be approximately written as Only downlink communications from the GBS and UAV to GUs are considered. The GBS and UAV are assumed to be equipped with an omnidirectional antenna, and they transmit data with constant power P G t and P U t . To facilitate the problem, we assume there exists only largescale fading for the GBS-GU links. The signal power that GU k receives from the GBS is modelled as where g k ≜ g 0 (d G k ) m accounts for path loss. g 0 denotes the reference channel power gain at d 0 = 1 m, while d G k = ‖q G − q k ‖ denotes the distance between the GBS and GU k, and m denotes the path loss exponent. u k is an independent and identically distributed log-normal random variable describing the shadow fading. Specifically, the decibel of u k follows a normal distribution with mean value 0, standard deviation , that is, 10 log 10 u k ∼ N (0, 2 ).
The UAV is assumed to communicate with GUs by LOS links, and Doppler effects caused by the UAV's motion is assumed to be well compensated. Therefore, the signal power GU k receives from the UAV can be described as where is the free-space path loss of the LOS channel, with h 0 denoting the channel power gain at the refer- − q k ‖ denotes the distance between UAV and GU k at time slot n.
The GBS and the UAV are assumed to utilise a common frequency spectrum with a bandwidth B, and the TDMA scheme is adopted to serve all GUs. Binary variables a G k [n] and a U k [n] are defined to describe the association between the UAV, GBS, and GUs. a G k [n] = 1 indicates that the GU k is served by the GBS in time slot n, and a U k [n] means similarly. We assume that the GBS and UAV can only serve one GU at a time separately. A GU cannot be associated with both the GBS and the UAV at the same time, too. Below are the constraints on user association schemes: The GUs are assumed to receive additive white Gaussian noise with power spectral density of N 0 . Taking co-channel interference into consideration, the signal-to-interference-plusnoise ratio (SINR) for GUs associated with the GBS is SINR for GUs served by the UAV is Therefore, we obtain the throughput of GU k in time slot n: a G k [n] and a U k [n] are coupled together in (17), making the user association problem difficult to tackle. Considering that there always exist users far apart from each other, the GBS and UAV would, respectively, choose these users to serve almost at all the time to achieve higher throughput. Therefore we make an assumption that Then a lower approximation of R k [n] is given as The corresponding lower approximation of the achievable average rate of GU k during the whole period T is

Problem formulation
The position of GBS and GUs, as well as channel state information, are assumed to be precisely known. For the sake of fairness, we choose the minimum average throughput among all users as our optimisation objective: Our goal is to maximise R min by jointly optimising the user association with UAV and GBS, along with the UAV's trajectory, under the above constraints. This problem is formulated as P1: Constraint (22b) is a more straightforward form of Equation (21). In the optimisation process, the maximum value of R min that is no greater than anyR L k is searched, which is in corresponds with the definition.

PROPOSED SOLUTION
(P1) is a mixed-integer non-convex optimisation problem that is hard to solve, so we relax the binary constraints on a G   [25] and partition the entire optimisation variables into two blocks, Then problem (P2) could be decomposed into two optimisation problems of user association and UAV trajectory. In the user association problem, we hold Q fixed and optimise A. In the UAV trajectory optimisation problem, we hold A fixed and optimise Q. Then, these two subproblems needed to be solved alternatively to approaching the optimisation solution if (P2).
However, the decomposed UAV trajectory problem is still challenging as it is a non-convex optimisation problem. Besides, R L k [n]'s logarithmic expression in (19) would consume much computation resources, leading to an increase in computing time. Therefore, we introduce a locally tight lower bound to substitute the objective and successively approximates the original objective with a sequence of convex approximations (SCA) [27].
We describe the whole optimisation algorithm with an efficient algorithmic framework called the BSUM and call it 'the relaxed BSUM algorithm'. The details are given in Subsection 3.1. As BSUM is proposed for large-scale data optimisation, our problem could be solved with high efficiency applying this framework. Furthermore, a parallel algorithm called parallel successive convex approximation is also included in the BSUM framework. Therefore, the BSUM makes it easier for us to change our current algorithm into a parallel version, also our next research target. Besides, our current algorithm faces the problem of local maximum points. Specifically, the algorithm may plunge into a certain point that is a local maximum point but not the global maximum point, and output a solution whose performance is much lower than the global optimisation value. But within the BSUM framework, we could consider more flexible variable update rules rather than what we use in the current algorithm. In fact, we are now in process of looking for a suitable update rule to jump out of local maximum points.
Apart from the method discussed above, we propose another way to deal with the local maximum points problem in Subsection 3.2. An improved solution that introduces a penalty factor and adjusts the iteration strategy is described in that subsection. In Subsection 3.3, the connections between the original problem (P1) and the transformed problems are given, and these two algorithms' convergence behaviours are discussed.

The relaxed BSUM algorithm
As has been stated, the optimisation problem (P2) is decomposed here into two subproblems: one to optimise the user association schemes assuming the UAV trajectory is known, and the other to optimise the UAV's trajectory with given user association schemes.

Subproblem of user association optimisation
With the binary constraint (14) relaxed and the UAV trajectory {q[n]} given, the user association question is as P3: This is a standard linear programming (LP) problem, which can be solved by today's mathematical solvers with high efficiency.

Subproblem of UAV trajectory optimisation
Assuming that the user association variables a G k [n], a U k [n] are given, the UAV's trajectory optimisation problem can be denoted as P4: With the definition of some constants and expressions, , (25b) can be rewritten in a more compact form: Constraint (26) where can be given by is a quadratic function of q[n] and always concave as l r k [n] is always negative. This is an attached benefit of our process of finding a lower concave bound of For ease of statements, a slack variable d k q [n] is introduced in place of the lower linear bound of ‖q[n] − q k ‖ 2 : As the two sides of (29) are always equal in the optimal point, d k q [n] is an equivalent substitution of the lower linear bound shown above.
Then the lower bound of E G k [n] is given as Till now, constraints (26) have been transformed into the standard convex constraints on arguments {q[n], d k q [n]}: However, in the existence of the logarithmic expression of E G k [n]'s lower bound, it is still high demanding to deal with the constraints (31). Therefore, we propose a method to approximate and calculate the logarithmic expression by replacing it with a concave fractional function: a(x 0 ), b(x 0 ), C (x 0 ) are coefficients calculated at local point x 0 , and we have carefully chosen their values to make g(x, x 0 ) as close to f (x) as possible. See the Appendix for details. Define , then ) .
The corresponding lower fractional bound can be expressed aŝ For simplicity, we denote a(x 0,k [n]) ⋅ ln 2, b(x 0,k [n]) ⋅ ln 2, C (x 0,k [n])∕ln 2 as a, b, C , and thenÊ G k [n] can be rewritten asÊ It could be verified with a simple analysis thatÊ G k [n] is concave concerning d k q [n]. Finally, the subproblem of UAV trajectory optimisation is approximated as a standard convex optimisation problem which can be efficiently solved by standard convex optimisation solvers, in the form of the following:

Initialisation algorithm
In our algorithms, the UAV's trajectory is optimised in an iteration way that needs an initial trajectory. Besides, the user association problem is also optimised, assuming a known trajectory of the UAV. Thus, an initialisation scheme of the UAV's trajectory is required as the beginning of our two algorithms. For simplification, we divide all GUs into two categories by power P G r,k of the signal they received from the GBS. Those GUs receiving higher P G r,k than P th are initially associated with the GBS while the others are linked to the UAV. Denote the indexes of UAV's users as where K U is the total number of UAV's users. Temporarily assuming that the UAV's altitude is fixed at H 0 , then the UAV and its associated GUs could be extracted as a simpler network: a UAV serving multiple fixed ground users. This simpler network has been discussed comprehensively in previous work like [29,30].
Intuitively, we consider two extreme cases concerning the UAV's flight period T : • (1) When T → 0, the UAV has to stay at one point q U 0 = (x 0 , y 0 , H 0 ) all the time. Best values of x 0 , y 0 could be obtained by solving the following problem: This problem aims to maximise the minimum communication rate R U min among all UAV's users by optimising x 0 , y 0 . Problem P6 is also non-convex, and we handle it with an iterative method that we have introduced in the trajectory optimisation subproblem. Limited by space, we omit the details to solve problem P6 and only show the initial point of iterations here.
• (2) When T → +∞, the UAV would have enough time that it can reach the top of each GU to enjoy the employ communication channel. As T is sufficiently large, we could neglect the UAV's time to fly among different users. Then T is divided between all users for the UAV to hover on. For fairness, the hovering time T i distributed to user k U i is inversely proportional to its communication rate R hov i : We obtain two trajectory initialisation schemes with low complexity based on the two extreme cases called 'circle-init' and 'TSP-init'. In circle-init, the UAV would fly around the point q U 0 , and the trajectory is designed as a circle centred on q U 0 . The circle's circumference is set as V h max T when T is small. As T increases, the circle's radius would reach r max , which is the maximum value for all circles centred on q U 0 and confined in the map.
In TSP-init, the UAV tries to fly over each associated GUs and hover on them. As the actual period T is limited, we have to minimise the UAV traveling time T trav to set aside more time for hovering. Considering that the UAV's maximum horizontal speed V h max T is constant, minimising T trav is equivalent to minimising the total travelling distance among different GUs. The latter is known as the classic travelling salesman problem (TSP), where a list of cities and the distances between each pair of cities are given, and the shortest possible route that visits

6:
until The fractional increase of the objective value is below a certain set threshold 1 .

7:
Round the continuous solution each city exactly once and returns to the origin city needs to be found. Although TSP is a non-deterministic polynomial complete problem, we could solve it efficiently by dynamic programming in our small-scale case (number of GUs connected to the UAV is not large). If T > T trav , we could take the TSP problem's solution as the UAV's flight trajectory. The hover time is divided like above: But if T ≤ T trav , then there is no time left for hovering. To design a closed trajectory, we scale down the TSP problem's solution to ensure that the UAV can exactly fly one round in the period T . The central point adopted as the reference point of re-scaling is q U 0 . As we derive circle-init and TSP-init from the two extreme cases, they are not universal. While circle-init does not suit the case when T is large, TSP-init does not suit when T is small. Therefore, we combine these two schemes as a new scheme named 'comb-init'. Comb-init is employed as the final initialisation scheme, where we invoke TSP-init when T > T trav , and invoke circle-init when T ⩽ T trav .

Overall algorithm
In the final, we obtain an iterative algorithm that successively solves the subproblems presented above by the BSUM method. For simplicity, we adopt the cyclic rule to select the variables to be optimised in each iteration. As problem (24) can be solved exactly, and problem (34) needs successively approximating, we first solve the user association problem for once and then conduct the UAV trajectory optimisation for k tra times in a cycle of iteration. Note that the solution we acquire by LP (24)

BSUM-based penalty algorithm
In this subsection, an algorithm applying the penalty factor based on the BSUM algorithm is proposed, called 'the BSUMbased penalty algorithm'. Meanwhile, a penalty factor is introduced to restrict a G k [n], a U k [n] to near 0 or 1.

Penalty factor
Unlike the first algorithm, a penalty factor is added to the objective function after the relaxation of the binary constraints. By introducing the penalty factor, a G k [n] and a U k [n] would be restricted to near 0 or 1 while the algorithm searches for more massive throughput. Therefore, this algorithm is more similar to the actual problem in the structure, where the values of aG[n], aU[n] are binary. Besides this, the changes in a G k [n], a U k [n] between different iterations also become much more extensive. In the relaxed BSUM algorithm, the searching areas of a G k [n], a U k [n] are mainly limited by the gradients of R min . When the penalty factor is applied, the objective function gradients are added with an item that inspires the algorithm to move towards those points whose a G k [n], a U k [n] are all 0 or 1, and helps the algorithm out of local maximum points. This advantage is validated in the numerical simulations. Function is chosen to form our penalty factors, with which we get a lower bound of R L k [n]: a lower bound of (40) is obtained as

Subproblem of UAV trajectory optimisation
In this subproblem, values of a G,r k [n], a U,r k [n] are all considered to be binary constants. Therefore, the subproblem of UAV trajectory optimisation with penalty factors is the same as the problem (25) and disposed of in the same method.

Overall algorithm
An iteration strategy similar to the BSUM algorithm is adopted here, which also alternatively solves user association and trajectory optimisation subproblems. A difference from the first algorithm is that the subproblem of user association is no longer a simple LP program, but added with the penalty factors and solved by the SCA approximation. Details are described in Algorithm 2.

15:
until The fractional increase of the objective value is below a certain set threshold 2 for 10 times.

Convergency analysis
This subsection analyses our algorithms' convergence behaviours and discusses the relationship between the original problem and the transformed problems. The first conversion is relaxation from the original problem (P1) to (P2). In this conversion, the objective function remains the same, represented as Γ (A, Q). Denoting the optimisation values of problems (P1) and (P2) as R ori , R rel , then there is R ori ⩽ R rel , as the relaxation operation enlarges A's and Q's feasible region. Then problem (P2) is partitioned into two subproblems (P3) and (P4) in both algorithms, whose optimisation values are lower bounds of R rel because the feasible area is narrowed down when A or Q is held.
In the relaxed BSUM algorithm, the user association subproblem (P3) is precisely processed by LP. If the value of Γ(A, Q) is Γ(A r , Q r ) before the rth iteration of user association optimisation, then it becomes Γ(A r+1 , Q r ) after this iteration. As this LP problem is globally solved, we have The UAV trajectory subproblem (P4) is approximated by the problem (P5), where the objective function is replaced bŷ Γ r (A, Q).Γ r (A, Q) is a global lower bound of Γ(A, Q), and satisfiesΓ Problem (P5) is also solved precisely, and the optimal objective value would therefore be no less than the function value at the Therefore, we have proved that Γ(A r , Q r ) ⩽ Γ(A r+1 , Q r+1 ). Because R rel is an upper bound of all subproblems' optimal objective value, we could know that Γ(A r , Q r ) would converge to a particular value R con as r increases. Influenced by the problem of local maximum points, there may be a massive difference between R con and R rel . Finally, the converged continuous solution is rounded to an integral solution A int , Q int . Meanwhile, the corresponding system throughput would fall to R int = Γ(A int , Q int ) ⩽ R ori .
In the BSUM-based penalty algorithm, there are several penalty-based optimisations after solving the LP problem in user association. These optimisations force variables A to be near binary, while they might lower the throughput, like the performance loss caused by the rounding operation in the relaxed BSUM algorithm. Therefore, this algorithm may keep vibrating rather than converge to one point. Fortunately, the simulation results indicate that this algorithm also roughly converge and that the vibration amplitude is pretty small.
Our analysis of the algorithms' convergency is validated in numerical simulations and visually described in Figure 2 with a concrete case. The optimal values of problems (P1) and (P2) are drawn schematically, as they cannot be obtained from our algorithms.

SIMULATION RESULTS
In this section, the system performance of our proposed schemes is evaluated via numerical simulations. A rectangular area of size 3 km × 3 km is considered, in which the GBS is fixed 20 m above the centre point. For comparison, the GUs and  channel parameters are generated by one random realisation of our program. All parameters used to construct the system and execute the algorithms are assigned as Table 1, except for those we specially stated in the following text or figures. All the simulations are conducted in the MATLAB environment applying the CVX, a package for specifying and solving convex programs [31].
To validate of our two proposed algorithms' effectiveness, they are compared with a fundamental scheme in which the user association scenarios of the GBS and UAV are optimised separately. In the fundamental scheme, GUs are divided into two groups, like in the initialisation scheme. The GUs with sufficient P G r,k are served by the GBS in turn, and the others are distributed to the UAV. Interference is considered in the fundamental scheme, while UAV's trajectory and user association scenarios are jointly optimised.

Iteration progress
The convergence behaviours of our proposed algorithms are described in Figure 2. The iteration times of trajectory optimisation that consumes the main part of computation time and resource is taken as the indicator to describe the algorithms' iter- The user association schemes ation progress. The blue curve depicts the continuous solution we obtain in the relaxed BSUM algorithm, while the red curve shows the rounded solution, that is, the final solution. Differences between these two curves reveal the rounding operation's performance loss, which is about 2%. The black curve denotes the BSUM-based penalty algorithm. We can conclude that the relaxed BSUM algorithm converges quickly and stably and that the BSUM-based penalty algorithm converges slowly and keeps fluctuating. However, the latter performs much better when it comes to max-min communication throughput. These results verified our discussion in Subsections 3.2 and 3.3. In the application, the relaxed BSUM algorithm is more appropriate for cases demanding low latency, while the BSUM-based penalty algorithm satisfies higher throughput requirements.

Optimisation results
In this subsection, the output results of different algorithms are compared in Figures 3 and 4. Intuitionally, the user association schemes and UAV trajectory generated by the relaxed BSUM algorithm and the BSUM-based penalty algorithm look alike. Therefore, we just display the results of the latter to save space. The results obtained from the fundamental scheme are also shown in comparison.
The user association schemes are represented in Figure 3, where the GUs are sorted in ascending order by the power  The unmanned aerial vehicle's (UAV) 3D trajectory P G r,k of signals they received from the GBS. In the fundamental scheme, the GUs are rigidly divided into two groups by P G r,k . The GBS chooses those GUs with higher P G r,k as its users, while the UAV serves the others. In this way, the UAV could assist those GUs who cannot communicate with the GBS efficiently. In addition, this method could also mitigate the co-channel interference as the GBS's interference on these GUs is quite weak.
In the BSUM-based penalty algorithm, the association schemes also adopt the division method, but more flexible. There are several GUs that receive service from both the GBS and the UAV. As a result of the joint optimisation, user switch behaviours of the GBS and UAV are related. For example, there is a short-term user switch for both of them at the time of near 40 s. This flexibility provides a considerable benefit for the improvement of system throughput, which would be shown later. Figure 4 depicts the UAV's optimised 3D trajectory along with its horizontal projection. All GUs are represented as coloured dots, which reveals the quality of their communication channels to the GBS. Those GUs with colour closer to green have relatively better channels to the GBS, while colour closer to red means the opposite. In a whole period, the UAV flies across the map to provide service for GUs in different positions. Roughly, the UAV tries to fly near those GUs with lower P G r,k while keeping a distance from the GU communicating with the GBS. When the UAV is close to its associated, it descends to utilise the best LOS channel. At other times, the UAV climbs up to avoid causing interference to the GBS's user.
The UAV's horizontal trajectories acquired from the BSUMbased penalty algorithm and the fundamental scheme are also compared. In Figure 4(a), the UAV's approaches the top of each associated GU in turn. However, in Figure 4(b), the UAV has to fly near the map border sometimes. The reason is that the GBS's user association schemes are not relevant to the UAV's trajectory or user arrangements in the fundamental scheme, and the UAV has to fly away if a nearby GU is communicating with the GBS. While in our BSUM-based penalty algorithm, the GBS's user association schemes are jointly optimised with other variables. Therefore, the GBS could select a GU far from the UAV for communicating, rather than keeping the UAV away.

Performance analysis
To validate our proposed algorithms' effectiveness, we first consider the fairness between different GUs via Figure 5. The average throughput of each GU is displayed as the green bars, and the minimum average throughput among all GUs is represented as the black dashed line. This figure is drawn using the BSUMbased penalty algorithm's outputs when the UAV's flight period T = 240 s, where the relative difference between the maximum and minimum average throughput values is 2.5%. When T is more prolonged, this indicator could descend to 1% due to the increased flexibility. The relaxed BSUM algorithm also shows a similar performance. These results demonstrate that our method of applying the minimum average throughput R min as the optimisation objective is quite beneficial to user fairness. Then the performance of our initialisation schemes circleinit, TSP-init, and comb-init are compared in Figure 6. As we predicted, circle-init outputs greater throughput when T is small, while TSP-init performs better when T is large. The combined algorithm comb-init accurately select and utilise the better algorithm in all tested cases. Besides, we also analysed the impacts of the UAV's height on the system performance by fixing the UAV's height in the algorithms. The results are displayed in Figure 7, where the 3D algorithms' performance is also shown as references. The height of the UAV influences the network in two aspects. On the one hand, a lower altitude could provide the UAV with better channels to communicate with GUs, thus is beneficial to the network. On the other hand, the users linked to the GBS would suffer more serious interference if the UAV flies down, which will worsen the system performance. These influences are apparent in the fundamental scheme, whose throughput first climbs up, but then falls down as the UAV's height increases. In our BSUM-based algorithm and the relaxed BSUM algorithm, the interference problem's influence is pretty reduced via the joint optimisation of trajectory and the UAV and GBS's user association schemes. Therefore, the impacts on the UAV's communi- cation channel would play a dominant role and keep improving the system performance as the UAV's height decreases to H min . Anyhow, the 3D schemes show better performance than 2D schemes in all algorithms. Figure 8 describes the throughput results obtained at different T . The performance of our two proposed algorithms and the fundamental scheme is contrasted with their 2D versions. Note that the UAV's height is also optimised to maximise the system throughput in the 2D schemes. In the BSUM-based penalty algorithm and the relaxed BSUM algorithm, the 3D versions only have small improvements upon the 2D versions. This is because the interference impact has been weakened a lot, as we stated above. Nevertheless, in the fundamental scheme, the 3D algorithm introduces a new optimisation dimension to relieve the co-channel interference, bringing a large gain. Moreover, this gain is enhanced as the increase of T , which allows the UAV to fly higher and more flexibly to avoid interference.
In the three 3D algorithms, the BSUM-based penalty algorithm behaves best, whose performance exceeds the fundamental one by around 20-30% in different cases. The throughput of the relaxed BSUM algorithm is similar to the fundamental scheme when T is small. As T increases, it improves rapidly and approaches the performance of the BSUM-based penalty algorithm.
Furthermore, the obtained max-min throughput grows with the increases of T in all algorithms, as the longer T provides more possible combinations of the user association schemes and a broader flight range. However, a longer T also means that a GU needs to wait more time to acquire service, and this is commonly referred to as the rate-delay balancing problem. When T is sufficiently large, the UAV could reach the map boundary area, and the user association scheme is near optimal, leading to the saturation of the throughput. Besides, the saturation values of the throughput in the 3D and 2D versions of the relaxed BSUM algorithms and the BSUM-based algorithms are close, while the 3D and 2D versions of the fundamental scheme have similar saturation values.

CONCLUSION
This paper proposes a design to reduce the interference in a UAV-assisted hybrid cellular network by choosing the suitable UAV trajectory and user association schemes. Specifically, the UAV and GBS's user association, along with UAV's 3D trajectory, are jointly optimised to maximise the minimum average throughput between all GUs with the existence of several constraints. By applying the framework of BSUM and introducing the penalty factor, two effective algorithms are, respectively, raised for application of low latency and QoS scenes. A method to simplify the calculation of logarithmic expressions in the optimisation problem is also put forward. Our proposed algorithms are numerically verified to be effective and outperform the fundamental scheme to a great extent.