Service-Aware User Association and Resource Allocation in Integrated Terrestrial and Non-Terrestrial Networks: A Genetic Algorithm Approach

In 6G networks and beyond, multiple radio access networks (RANs), including; the satellite, high altitude platforms, low altitude platforms, and the terrestrial network, will co-exist. These networks are characterized by different capabilities and limitations in meeting the envisioned 6G contrasting user requirements. Therefore, associating users with the appropriate radio access network (RAN) in such an integrated network is rigorous and complex. In this work, the user association and resource allocation problem is formulated as a multi-objective optimization problem (MOOP), aiming to maximize data rate while minimizing mobility-induced handoff in the integrated network. Moreover, the problem is formulated in such a way as to prioritize the service provisioning of mission-critical users. The weighted sum method is adopted to simplify and transform the MOOP into a single-objective optimization problem (SOOP). In order to solve the formulated NP-hard SOOP, a genetic algorithm (GA) whose fitness value is based on the user’s service group is proposed. The performance of the proposed algorithm is evaluated by comparing it to the optimal solution, the greedy signal-to-interference-plus-noise ratio (SINR) based association, and the random user association algorithms. Simulation results show that as the number of access nodes in the network increases, the GA’s spectrum efficiency (SE) remains within 0.4% of the optimal solution. Moreover, the GA outperforms all three schemes in user acceptance ratio (AR) and handoff probability.


I. INTRODUCTION A. BACKGROUND
In recent times, the information communication technology sector has witnessed an explosive growth in demand for high-speed wireless access, which has increasingly strained the terrestrial network (TN) [1]. While technologies such as ultra-dense networks (UDNs) and device-to-device (D2D) communications have shown great potential in increasing the The associate editor coordinating the review of this manuscript and approving it for publication was Xiaodong Xu . capacity of the terrestrial networks (TNs), they are faced with different challenges [2], [3]. UDNs are limited by frequent handoff, interference, and backhaul challenges, while D2D communication is faced with frequency planning and resource allocation implementation complexity.
To solve such challenges and increase the TN's capabilities in providing ubiquitous broadband connectivity, one of the key enabling features of the sixth generation of wireless networks (6G) is the integration of TNs with non-terrestrial networks (NTNs) [4], [5], [6], [7], [8]. The considered NTNs include satellite communications (SatComs) and unmanned VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ aerial vehicles (UAVs), which can be classified into high altitude platforms (HAPs) and low altitude platforms (LAPs) [2]. Fig. 1 depicts the integrated terrestrial and non-terrestrial network (ITNTN) architecture, which is a layered and 3-dimensional integration of SatComs, aerial platforms, i.e., HAPs and LAPs, and the TN. SatComs consists of the Geostationary Earth Orbit (GEO), Medium Earth Orbit (MEO), and Low Earth Orbit (LEO) located at altitudes of 35786 km, 7000 -20,000 km, and 600 -1500 km, respectively [9]. On the other hand, HAPs are repeaters flying at an altitude of 17-22 km in the stratosphere [10]. They can be classified as aerostatic and aerodynamic, with aerostatic HAPs taking the form of either balloons or airships. Balloons are designed to stay still in space, while airships are quasi-stationary with onboard electric motors and propellers for station keeping [10]. On the contrary, the aerodynamic HAP is an aircraft that has to stay in a forward motion to keep afloat. Compared to the other networks of the ITNTN, LAPs have the lowest cost and are characterized by fast and easy deployment [10], [11]. This attribute makes them suitable for providing communication services for emergency response and acting as aerial BSs for direct user equipment (UE) connectivity and traffic offloading during limited-duration events such as festivals and sports events. They are relatively small and light and operate at low altitudes, not exceeding 10 km above the earth's surface. The TN layer includes fixed and mobile users using radio access technologies such as 5G, 4G, and WiFi to access various heterogeneous networks consisting of macro, micro, pico, and femtocells. The integration of TN with NTNs is motivated by the technological advancement in manufacturing and launching processes that has resulted in the massive deployment of LEO and MEO satellites by companies such as OneWeb, and SpaceX [12]. Furthermore, companies such as Airbus and Google have heavily invested in HAPs [13], [14]. Moreover, the third generation partnership project (3GPP) has successfully conducted feasibility studies on radio access through the NTNs [9], [15]. According to 3GPP, the NTNs will support the TN in providing radio access via two radio access network (RAN) architectures depicted by Fig. 2. First is the transparent/bent pipe payload architecture portrayed by Fig. 2(a), in which the NTN access node (AN) acts only as a relay, with its function being frequency filtering, frequency conversion from uplink to downlink, and signal amplification. The second architecture is the regenerative payload configuration illustrated in Fig. 2(b) in which the AN not only performs the functions of a transparent payload but also implements demodulation/decoding, switch and/or routing, coding/modulation. Such an AN is considered a next-generation node B (gNB) as it performs the functions of a base station. In both architectures, the UE can connect to the NTN AN either directly or through a relay node (RN) using the new radio (NR) interface, as depicted by Fig. 2. On the other hand, the connection from the transponder to the gateway (GW) is through the NR air interface for the transparent configuration, while the regenerative architecture uses the NG (Next Generation) interface. An interested reader is referred to [9] for a detailed description of the NTN RAN architecture. Indeed, the future wireless networks will have the TN co-existing with the NTNs to provide radio access to multi-mode UE.
Motivated by these trends, this paper seeks to address the problem of user association and resource allocation in the ITNTN. In particular, and similar to [16] and [17], the work considers the usage scenario in which there exists a large number of users whose traffic cannot be entirely supported by the TN RAN, for example, in an urban area during a carnival event. Such a usage scenario necessitates the deployment of NTNs to de-congest and support the TN RAN in providing radio access to the different users. The UE in this integrated network is considered a multi-radio terminal with the ability to access either the terrestrial access network or the space/airborne communication networks, including the SatComs, HAPs, and LAPs. Without loss of generality, the NTN access nodes (ANs) are assumed to be regenerative, and users can connect to them directly. The different radio access networks (RANs) in this ITNTN have diverse constraints and capabilities towards meeting the heterogeneous requirements of the beyond 5G (B5G) users.
One of the salient features of the fifth generation of wireless networks (5G) is differentiated service provisioning. Unlike the previous generations of mobile networks that were based on the traditional one-size-fits-all architecture, 5G leverages technologies such as network function virtualization (NFV) [18], software defined networking (SDN) [3], and network slicing [19] to provide an agile and flexible network. Such a network is designed for service provisioning based on the heterogeneous quality of service (QoS) requirements of diverse users. To this end, the various user services are classified into different service groups (use-cases), with each service group comprising services of related attributes and priorities [20], [21]. The envisaged service groups in 6G according to [5] include further enhanced mobile broadband (feMBB), enhanced ultra-reliable and low-latency communications (euRLLC), ultra-massive machine-type communications (umMTC), long-distance and high-mobility communication (LDHMC), and extremely low power communications (ELPC).
The problem that arises then is how to map such heterogeneous requirements to the appropriate RAN in the ITNTN. For example, although the TN is endowed with a rich pool of resources that result in high throughput and low latency, it still cannot support an increased number of users in situations of traffic overload; hence it is essential to decide which users are served by the TN and which by the NTN RANs. Moreover, the TN is limited in coverage per base station (BS), making it an unsuitable RAN for LDHMC users. In comparison, the NTNs are characterized by a wide coverage that reduces the number of handoffs experienced by the LDHMC users, ultimately decreasing the associated delays and signaling overheads. Along the same lines, SatComs is limited by a high end-to-end delay and is not suitable for direct connectivity of the mission-critical euRLLC service group. Similarly, the HAPs and LAPs are constrained by the number of available channels; hence, it is crucial to prioritize their use for traffic whose QoS requirements may not be met satisfactorily by the other RANs. Clearly, associating the different users to appropriate RANs in the ITNTN while maximizing resource utilization and providing the required user QoS is a rigorous and complex task.

B. RELATED LITERATURE AND CONTRIBUTIONS
Several works have addressed resource allocation in the ITNTN in recent times. The authors in [22] propose a user association and resource allocation problem that maximizes the total data rate of an integrated UAVs and terrestrial RAN together with satellite and macro BSs backhaul links. This work is supplemented in [16] by the same authors by taking into account the access node's energy cost and the limited life endurance of UAVs. However, the work does not consider the satellite as an AN but only as a backhaul link and does not consider HAPs in the integration. Therefore, the authors did not consider the idea of having space (satellite) and aerial (LAPs and HAPs) nodes complement the TN in providing radio access to users. Moreover, the work did not account for how the different user requirements were met.
In [23], the authors propose a load balancing algorithm for an integrated satellite and TN. They define the radio resource utilization ratio as a metric used to measure each cell's load status and group traffic into delay-sensitive and delay-tolerant. Traffic is offloaded from an overloaded terrestrial cell to neighboring terrestrial cells first and then to the satellite cell. Delay-intolerant traffic is not offloaded to the satellite network. However, for 6G networks to have an effective and efficient resource allocation scheme, traffic can not only be classified into two but in different use-cases that effectively capture all future network demands. Moreover, this work only associates users to the NTN when the TN is overloaded. Besides, the authors did not consider other NTNs such as LAPs and HAPs.
In [24], the authors jointly optimize user association and resource allocation of both access and backhaul links, together with HAPs' locations, to maximize users' throughput in an integrated satellite, airborne, and TN. These authors neither considered the heterogeneous user QoS requirements nor the different limitations of the networks in the integration. The authors in [25] propose a joint algorithm that optimizes user association and resource allocation to a terrestrial macro base station (MBS) and multiple UAVs mounted BSs using in-band wireless backhaul. However, the authors only consider UAVs and TNs in their analysis, and they do not account for the provisioning of the different service group QoS requirements.
The authors in [26] propose an algorithm that maximizes the energy efficiency (EE) of an integrated satellite/terrestrial cache-enabled RAN. Both terrestrial ANs and LEO satellites provide content distribution and retrieval services, thereby offloading such traffic from the MBS. The work does not consider the implementation of an integrated terrestrial-aerialspace RAN. It neither addresses QoS provisioning to different use-cases, as it considers only one service: content distribution. In [17], the authors propose an optimization problem that maximizes the network data rate while ensuring QoS by minimizing interference in an integrated satellite and TN. The work fails to distinguish users according to their different QoS requirements, such as data rates, and does not guarantee that these QoS demands are satisfied. Besides, it assumes that all the available networks can support all different users.
The authors in [27] maximized the minimum ergodic achievable rate of a user-UAV link. This work optimizes only the UAV RAN and not the entire integrated terrestrialsatellite-UAV RAN. Besides, the work did not incorporate QoS provisioning of the different use-cases. The authors in [28] first optimize multi-beam dynamic radio resource allocation for LEO-ground downlinks. Next, they optimize the dynamic resource allocation for HAP-ground downlinks when HAPs and LEO satellites share the same spectrum. This work does not consider joint resource management of the different RANs, as it optimizes the resource allocation of each RAN separately. The authors in [29] propose a traffic offloading scheme that considers the diversity of user demands. The optimization algorithm maximizes the eMBB users' data rate, subjected to stringent outage probability of the uRLLC usecase. This work does not consider using SatComs directly for radio access, and neither does it incorporate the UAVs (i.e., LAPs and HAPs) in the radio resource management.
Different from all the above work, we seek to find an optimal dynamic radio access user association and resource allocation scheme for an ITNTN comprising the TN, LAP, HAP, and SatComs ANs. We consider the heterogeneity in users' demands and the different RANs' uniqueness in meeting these diverse demands. The work considers three service groups: feMBB, euRLLC, and LDHMC.
The user association and resource allocation problem is formulated as a multi-objective optimization problem (MOOP), maximizing the total network data rate while simultaneously prioritizing large coverage NTNs over the TNs for service provisioning of the mobile LDHMC service group, with the objective of minimizing mobility-induced handoff probability. Moreover, given that denial of service to the euRLLC service group may result in catastrophic events, the optimization problem prioritizes the service provisioning of this service group over others and limits its access to the SatComs AN. The MOOP is simplified and transformed into a weighted sum single-objective optimization problem (SOOP). The formulated combinatorial and non-convex SOOP is NP-hard, with no efficient polynomialtime solution. Owing to its simplicity and efficiency in solving non-convex and combinatorial problems [30], we solve the SOOP using the genetic algorithm (GA).
The GA is a meta-heuristic search algorithm that uses the theory of evolution and natural selection to solve optimization problems. It is well suited for multi-objective and nonmathematical optimization problems, efficiently and easily enforcing different constraints, and also searching over multiple sets of solutions in a large search space to return a near-optimal solution [31]. Given its ease of implementation and optimization of discrete and continuous radio parameters, the GA is an excellent optimization tool for radio resource management in the ITNTN. The performance of the proposed GA is compared to three other algorithms; the optimal solution simulated using the gurobi solver, the greedy signalto-interference-plus-noise ratio (SINR) based solution [32], and the random user association (RUA) solution. Simulation results show that as the number of ANs increase in the ITNTN, the GA's spectrum efficiency (SE) performance remains within 0.4% of the optimal solution and outperforms the greedy and RUA by 1.23% and 0.97%, respectively. Moreover, the GA outperforms the optimal, greedy, and RUA algorithms in handoff probability by 8.4%, 14.9%, and 51.8% on average, respectively. Furthermore, the GA shows an acceptance ratio (AR) performance that is better than the optimal, the greedy, and the RUA solutions by 1.41%, 10.8%, and 7.6%, respectively. The key contributions of this work can therefore be summarised as follows: • Formulation of a user association and resource allocation optimization problem that maximizes the data rate of the ITNTN while simultaneously minimizing the probability of mobility-induced handoff. Handoff is minimized by prioritizing the use of large coverage NTNs by the mobile LDHMC service group. Moreover, service provisioning of the mission-critical euRLLC use-case is prioritized over other use-cases.
• The formulated multi-objective problem is transformed into a single-objective problem which is solved using the GA by encoding the problem into a sequence of chromosomes, with the genes representing user association solutions. Service group-dependent fitness functions are formulated to determine the near-optimal user association and resource allocation solution.
• Numerical results are presented, comparing the proposed GA to three algorithms; the optimal solution simulated using the gurobi solver, the greedy algorithm whose association is based on maximum SINR [32], and the random user association solution.

C. ORGANISATION
The remainder of this paper is organized as follows: Section II gives a detailed description of the system model and assumptions, while section III defines the problem formulation. In Section IV, the GA process is reviewed, and the solution to the formulated problem based on the GA described. The results are analyzed in section V, and finally, the conclusion presented in section VI. The notations used in this paper are presented in Table 1.

II. SYSTEM MODEL
This section describes the deployment scenario, mobility model, channel model, signal quality model and assumptions considered in this work.

A. DEPLOYMENT SCENARIO
A downlink transmission of an integrated communication network consisting of four RANs namely, the MBS, LAP, HAP and the LEO SatComs as depicted in Fig. 3 [33], [34], [35], and [36], this work is premised on the assumption that each UAV AN o j ∈ {u, h} is stationary or quasi-stationary with negligible mobility. Immobility of UAVs is assumed to avoid disconnections due to the UAV AN moving out of coverage of already connected users, some of whom could be mission critical. Consequently, we assume to use the rotary-wing LAPs and balloon or airship HAPs that have the ability to remain quasi-stationary [11], [37]. Moreover, we assume that the placement of the UAV ANs has already been optimized to cater to the usage scenario in which there is a large number of users, say in an urban area, during a carnival event.
The system provides downlink communications to a set of users I = {1, 2, 3, . . . , |I|} and supports three use-cases namely; feMBB, euRLLC, and LDHMC, denoted by E, R, and D respectively, such that υ ∈ {E, R, D}. Users demanding use-cases E, R, and D are grouped in sets denoted by I E , I R , and I D respectively, such that, i E ∈ I E , i R ∈ I R , i D ∈ I D , I = I E ∪ I R ∪ I D and I E ∩ I R ∩ I D = ∅. Without loss of generality, the I E and I R users are assumed to be static while the I D users are mobile. An example of a static use-case belonging to the I R service group is remote surgery within a hospital or private clinic in the considered urban area. Also, a user i ∈ I is assumed to be embedded with multiple RAN interfaces, and thus can access any of the available RAN j ∈ {B, U, H, S} within its coverage.
Since the different RANs in the ITNTN may have different multiple access schemes such as OFDMA, TDMA, and FDMA, the basic bandwidth unit (BBU) is used to represent the unit of radio resources as was done in [38]. Therefore, VOLUME 10, 2022 no matter what type of access technique is used, the system capacity is represented in terms of bandwidth. Similar to [39], this work assumes shared spectrum mode, such that all ANs that belong to the same RAN j own the same set C j of BBUs.
Furthermore, as was done by [24], we make the following assumptions: i) A user i ∈ I can associate with at most one AN. ii) For simplicity, the different RANs use frequencies sparsely separated from each other, thus have no crosstier interference. iii) Intra-cell interference on the downlink between users associated with the same AN is negligible, as it can be effectively controlled through multiple access techniques. iv) Co-tier interference exists, where a user receives signals from different ANs in the same RAN. This interference is added to the thermal noise in the SINR expression defined in (9).
The system model used for the LEO SatComs RAN is adopted from [26]. A LEO satellite follows a specific pattern in which it periodically serves one area followed by another. Consequently, in this work, we assume that each non-overlapping area depicted in Fig. 3 is served by one LEO satellite for a given duration of time. Hence, a user in a given location can access only one satellite at any given time. As an illustration, let us consider each LEO satellite's view time of any given area to be t. Then if LEO 2 starts to view area 2 at time t1, this view continues until a time t1 + t, after which service provisioning of area 2 is handed over to LEO 1, which views from t1 + t until t1 + 2t. Handover from one satellite to another is assumed to be managed by the Network Control Center [26]. In this work, this handover process is out of scope and thus shall not be considered. As the entire considered terrestrial area is under coverage by a single satellite during its service period, this implies that seamless coverage of the terrestrial RAN is guaranteed by the satellite, irrespective of its mobility. Moreover, in practice, thousands of LEO satellites are deployed, for example, in the Starlink project, implying that a particular area can be covered by multiple LEO satellites simultaneously [40], [41]. Therefore, the mobility model of the satellite shall not be considered in this work.

B. USER MOBILITY MODEL
Due to its simplicity and analytical tractability, the random walk mobility model [42] is used to imitate the LDHMC users' movement patterns. In each new interval, a user i D ∈ I D chooses a direction θ i D ∈ [0 2π] that is randomly and uniformly distributed. In the same manner, the user's speed is randomly assigned following a uniform distribution, with V min and V max being the minimum and maximum velocity respectively, that a user can have. Considering the user's location in time interval t as ( where D max is the maximum distance that can be moved by a user in a given flight interval. We consider LDHMC users to move within the LEO satellite coverage area, assumed to be 5 km, such that they are reflected off the boundary of the circular region.

C. CHANNEL MODEL
The channel modelling is similar to our work in [32]. The ground location of any AN o j ∈ {B∪U ∪H∪S} is represented by In this case, for the NTN ANs, t o j is their projection on the ground. On the other hand, the coordinate of a user i ∈ I is given by where z o j is the height of the AN and . 2 is the 2-norm operator. Path loss modeling is divided into three categories: (i) MBS TN, (ii) UAV, i.e., HAPs and LAPs, and (iii) SatComs.

1) MBS TERRESTRIAL NETWORK PATH LOSS MODEL
The Path loss of a user i ∈ I using a BBU c j ∈ C j |j ∈ B to communicate to an urban MBS o j ∈ B that is located a distance d o j ,i meters away is given by [43] The term h o j represents the MBS height in meters, f o j ,c j is the carrier frequency in MHz, and τ o j ,i is the path loss due to shadow fading, assumed to be a Gaussian random variable with zero mean and σ standard deviation in dB. τ o j ,i can be expressed as τ o j ,i = log 10 (F o j ,i ) where F o j ,i is the log-normal shadow fading path loss between the user i and AN o j [43].

2) UAV PATH LOSS MODEL
The path loss from a UAV AN o j ∈ {U ∪ H} to a user i ∈ I is modelled according to [44]. In this model, the probability that a user i has a line of sight (LoS) link from a UAV AN o j ∈ {U ∪ H} is given by The constants x and y are dependent on the environment while the elevation angle θ o j ,i is given by . The signal from the UAV AN propagates through free space until it reaches the environment on earth, where it undergoes additional loss due to shadowing, scattering, and reflections caused by buildings, foliage, etcetera. Consequently, the path loss from the UAV AN o j ∈ {U ∪ H} and a user i located at distance d o j ,i is given by where d o j ,i is distance in meters computed using (2), f o j ,c j is carrier frequency in MHz, and η LoS or η NLoS is the additional loss experienced in LoS or NLOS propagation respectively. The addition loss η LoS or η NLoS has a Gaussian distribution [44]. The equivalent path loss is then given by [25] PL equiv where Prob NLoS o j ,i = 1 − Prob LoS o j ,i is the probability that a user experiences a NLoS link.

3) SATELLITE COMMUNICATIONS PATH LOSS MODEL
Assuming clear sky conditions, the path loss between an AN o j ∈ {S} and a user i at a distance d o j ,i from the node, is given by [9] PL o j ,i,c j = 20 log 10 d o j ,i + 20 log 10 The frequency f o j ,c j is in MHz, while the range d o j ,i in meters and can be determined using (2). On the other hand, τ o j in (7) depicts the loss due to shadow fading, while CL gives the clutter loss resulting from reflections and scattering caused by surrounding buildings and objects on the ground. τ o j and CL are dependent on whether the propagation is LoS or NLoS. The terms PL k , PL e , and PL y represent attenuation due to atmospheric gases, ionospheric or tropospheric scintillation, and building entry loss. The equivalent path loss is determined using (6) where by the LoS probability P LoS o j ,i which depends on the elevation angle and UE environment is obtained from Table 6.6.1-1 in [9].

4) SIGNAL QUALITY MODEL
Using the results from the preceding sections II-C1 to II-C3, the linear channel gain due to large scale fading effects of path loss and shadowing from an AN o j ∈ O and user i, communicating via a BBU c j is given by [43] In practice, for time division duplex systems, the ANs can exploit reciprocity between downlink and uplink channels to estimate the downlink channel [45], [46]. In comparison, for the frequency division duplex systems, there often exists weaker reciprocity in the uplink and downlink frequencies [45], [46]. Consequently, the channel state information (CSI) at the ANs for these systems can be obtained as feedback from the UE [45], [46]. Several works in literature are dedicated to reducing the amount of CSI feedback in the FDD systems [46], [47], but this is out of the scope of this work. The obtained CSI feedback from UEs suffers high latency for ANs at very high altitudes, like the SatComs; hence, this CSI is usually outdated [48]. Therefore, for such RANs, CSI can be obtained by utilizing the widely used training data-based CSI estimation techniques [49], [50], [51]. However, due to the complexities involved in simulating the CSI feedback and estimation techniques, this work assumes that the large-scale CSI is available at the ANs as was done in [24], [43], [52], and [27]. Therefore, the corresponding channel gain can be calculated using (8).
The SINR ratio experienced by the user i using BBU c j is then determined using The terms β o j ,iq,c j and β l j ,iq,c j denote the small-scale fast fading component that accounts for mobility of a user and is assumed to be independent and identically distributed (i.i.d) as CN (0, 1) [53]. P o j ,iq,c j is the transmit power from an AN o j to a user i q using BBU c j . Similar to [54], [55], and [56], VOLUME 10, 2022 and in a bid to reduce complexity, this work shall not optimize power allocation. Therefore, the transmit power P o j ,iq,c j is assumed to be fixed and uniformly allocated to all AN's BBUs. The power P o j ,iq,c j is given as P thres o j /|C j |, where P thres o j is the maximum available power at the AN o j and |.| represents the cardinality of the set. P l j ,iq,c j depicts the co-tier interference from any other AN l j in the j th RAN reaching the same user i q .
The maximum data rate that can be achieved by a user i ∈ I on the BBU c j of AN o j is then given by the Shannon capacity as where T c j is the bandwidth for BBU c j ∈ C j . Ultimately, the overall system data rate is then expressed by where ρ iυ ∈ [0, 1] is a weighting factor for a user i ∈ I demanding a given service belonging to service group υ ∈ {E, R, D}. ρ iυ prioritizes service provisioning of users based on the service they demand. µ o j ,i ∈ {0, 1} is a binary user association variable that specifies whether the user i ∈ I is associated with the AN o j ∈ O of the j th RAN, and is defined as: On the other hand, ω o j ,i,c j ∈ {0, 1} is a binary resource allocation factor that is 1 if the BBU c j ∈ C j of AN o j is allocated to the user i, and 0 otherwise, that is; We assume that an AN's BBU can be allocated to at most one user, as illustrated by (14).
Also, a user i ∈ I can only associate with ANs in whose coverage the user lies. Therefore, we define a binary index π o j ,i that indicates whether a user i is within AN o j s coverage or not. If the distance between a user i s location (x i , y i ) and the ground location of the AN o j (x o j , y o j ) is less than the AN's cell radius R o j , then π o j ,i = 1, otherwise, π o j ,i = 0, as illustrated by (15).
This work maximizes the overall system data rate while simultaneously minimizing the probability of mobilityinduced handoff. We define the probability of handoff as the ratio of the number of users that experienced a handoff to the total number of mobile users during a given transmission time interval (TTI). To minimize the probability of handoff, a handoff reduction function given in (16) is maximized in each TTI.
The term ζ in (16) represents the largest cell radius in the ITNTN. A mobile user traverses several small cells during an ongoing call. Hence, the smaller the ANs' cell radii, the more the number of handoffs due to user mobility. In order to limit the mobility-induced handoff, we introduce the ratio (R o j /ζ ) ∈ (0, 1] in (16) to prioritize the association of mobile users to ANs with the largest cell radius. In this way, a maximum value of the ratio (R o j /ζ ) implies a reduced number of handoffs since the mobile user will be associated with an available AN having the largest cell radius in each TTI, ultimately reducing the probability of handoff.

III. PROBLEM FORMULATION
In this section, the user association and resource allocation problem is formulated as a MOOP that maximizes the total system data rate while at the same time minimizing the probability of mobility-induced handoff. The probability of handoff is minimized by maximizing the association of mobile users to ANs with a large cell radius. Moreover, ρ iυ in (11) and (16) can be used to prioritize the service group R users over other users. The MOOP is formulated as where is the set consisting of all feasible user association and resource allocation solutions satisfying the following constraints Constraint C1 ensures that a user can only associate with an AN in whose coverage radius the user lies. C2 indicates that an i R user can be served by only one of either the MBS, LAP, or HAP AN at a time. C3 ensures that an euRLLc user is not attached to the satellite AN. In constraint C4, a feMBB or LDHSC user can associate with only one of the available ANs at any given time. Constraint C5 guarantees served users a minimum QoS in terms of data rate. I served is a set of users for which π o j ,i µ o j ,i ω o j ,i,c j = 1 while L υ thres is the minimum required data rate of a user demanding service group υ. In C6, the total radio resources allocated to all users attached to an AN o j must not exceed the AN's radio resource budget o j . Constraint C7 ensures that an AN's BBU can only be allocated to at most one user, while C8 gives the decision variables that are binary in nature.
To solve the MOOP in (17), the concept of Pareto optimality as defined in [57] is utilized, and this states that: Definition 1: A point 0 ∈ is said to be Pareto optimal if and only if there is no other point 1 ∈ such that ( 1 ) ≥ ( 0 ), ( 1 ) ≥ ( 0 ) and at least one or has been strictly improved.
In simple terms, a point is Pareto Optimal if there is no other point that can improve both and simultaneously. The set of all Pareto optimal points gives an optimal trade-off between and by providing the maximum value of for any given value of and vice verse. The weighted sum method is capable of providing a complete set of Pareto optimal solutions to the MOOP in (17). Therefore, similar to [58], [59], and [60] and owing to its simplicity and low computational complexity; the weighted sum method is adopted to transform the MOOP in (17) into a SOOP, which is a weighted sum of the two objective functions as shown below: 17a, 17b, 17c, 17d, 17e, 17f , 17g, 17h.
The terms δ 1 and δ 2 in (18) are normalization factors associated with and respectively, introduced to put the two functions on the same scale. α in (18) is used to provide a trade-off between data rate maximization and mobility-induced handoff minimization for the mobile users. This parameter is particularly useful since the ANs that maximize data rate may not necessarily minimize handoff, and hence is used to define the importance of the two objective functions. It is important to note that for static users, i.e., the euRLLC and feMBB service groups, α is 1 since these users do not experience mobility-induced handoff. On the other hand, when the values of α are varied in the range [0,1] for the mobile users, the set of all Pareto optimal points to the MOOP in (17) can be obtained [61]. However, in this work, we consider the priority of the mobile LDHMC service group to be handoff minimization, and as such, we set α to 0 for these users.
The optimization problem described in (18) is non-convex and combinatorial, making it NP-hard with no efficient polynomial-time solution. Exact algorithms such as branch and bound exist that return a global optimal solution for such problems. However, these algorithms have to search for all possible solutions (in the worst-case scenario) in the search space. Consequently, the exact algorithms have a high computational complexity that increases exponentially with the network's number of ANs and users. However, the GA has been proven to give near-optimal solutions with a polynomial-time complexity to non-convex and combinatorial problems [30]. Consequently, in this work, we adopt the GA to solve the problem formulated in (18).

IV. PROPOSED SOLUTION
The proposed GA solution is described in this section, but first, a brief review of GA is elucidated.

A. A REVIEW OF THE GENETIC ALGORITHM (GA)
The GA is a search meta-heuristic based on the principle of natural selection in which the fittest individuals of a population are selected to produce children of the next generation. The algorithm starts with generating an initial population consisting of a set of randomly generated solutions, also called chromosomes. Each chromosome is made up of genes, which are essentially the decision variables of the optimization problem.
A fitness function corresponding to the objective function is defined and used to measure the fitness of each chromosome in the population. Parents are selected from the population for reproduction based on their fitness values. In the crossover phase, the parents exchange genes, producing new chromosomes, otherwise called children. The crossover phase is governed by the probability of crossover P c , which determines whether to consider a child or parent chromosome in the new population.
The mutation operator is used to randomly change one or more genes of the chromosomes to create diversity in the new population and prevent the GA from converging to a local optimum. The mutation is dependent on the probability of mutation P m , with high values of P m changing the algorithm to random search. After mutation, the elitism operator ensures that the best solutions/chromosomes in the old population are not lost through crossover and mutation. Therefore, depending on the elitism ratio E r , a given fraction of the best chromosomes in the old population replace other chromosomes in the new population. The algorithm then terminates when the termination criteria are satisfied, and the chromosome with the best fitness value is the solution to the optimization problem. There are several termination criteria used in literature [62], including; (i) when a fixed number of generations is reached, (ii) when a certain fitness level is reached, and (iii) when there is no improvement in the best fitness value. In this work, the termination criterion is either when there is no improvement in the best fitness value for a given number of consecutive iterations or when a fixed number of generations is reached. The GA process is illustrated in Fig. 4.

B. THE PROPOSED GA
A population set M consisting of M chromosomes is defined. A chromosome µ k ∈ M is defined as an |I| dimensional vector that represents a user association solution In any given iteration, the fitness value of a gene µ k o j ,i of a chromosome µ k , that represents the optimization problem in (18), is defined depending on the service group of user i.

1) FITNESS VALUE FOR AN euRLLC OR feMBB GENE (DATA RATE MAXIMIZATION FITNESS VALUE)
The fitness value of a gene belonging to an i R or i E user is given by where δL o j ,i,c j is the normalized data rate a user i can achieve over one BBU c j of AN o j . x i,k ∈ {0, 1} denotes the validity status of the gene µ k o j ,i . x i,k = 1 if the user i s association variable/gene µ k o j ,i is valid and 0 otherwise. A gene µ k o j ,i is valid if the AN o j is within coverage of user i, is capable of serving the user, and has sufficient resources to meet the QoS requirements of the user. ξ i υ ∈ [0, 1] is the penalty cost for not admitting a user of service group υ ∈ {E, R}.

2) FITNESS VALUE FOR AN LDHMC GENE (HANDOFF MINIMIZATION FITNESS VALUE)
The fitness value for users demanding service group D prioritizes association to ANs with large cell radius in order to minimize handoff probability. Consequently, the fitness value for these users is expressed as where ξ i D is a penalty for not admitting an i D user. The penalties ξ i υ in (19) and ξ i D in (20) are varied depending on the priority of a given service group. We prioritize the euRLLC service group in this work since denial of its service may lead to catastrophic consequences. We also prioritize access to the NTN with a large coverage area for the mobile LDHMC use-case to minimize the probability of handoff. Nonetheless, the priority can also be based on other factors, such as the use-case that yields more revenue to the operator.
For any given gene µ k o j ,i of a chromosome µ k , if x i,k = 1, then the user i is allocated its required number of BBUs N BBU i given by (21), else, the user is not allocated any resources. We assume that L υ thres is the data rate request from a user demanding a service in the service group υ ∈ {R, E, D}.
The overall fitness value of the chromosome µ k in any given iteration is the summation of all fitness values for the different genes in the chromosome given by: According to (22), all users in the network cooperate and contribute to the fitness value of the chromosome, ultimately leading to fairness between users since the goal is to admit as many users as possible to maximize the fitness value.
The pseudo-code of the proposed GA is given in Algorithm 1. The selection of parents in line 7 is achieved using the roulette wheel technique. In this technique, each chromosome in the population is assigned a probability P µ k of being selected depicted by (23) that is proportional to its fitness value. The chromosomes with higher values of P µ k have higher chances of contributing to the creation of the next generation.
In lines 8-9, the two selected parents recombine through the crossover operator. In this work, the two-point crossover technique is used [62]. In this technique, two points are randomly selected on both parents, and genes are exchanged between them to create two different chromosomes, otherwise called children. For each created child, if a randomly generated number in the range [0, 1] is less than P c , then the child is inserted into the new population; otherwise, the parent is.
The new chromosomes are then mutated as per line 10 of Algorithm 1. For each gene in the chromosomes created in lines 8-9, if a random number in the range [0,1] is less than P m , the gene is mutated by replacing it with a random gene µ k o j ,i ∈ {1, 2, . . . , |O|}, otherwise, the gene is not mutated. Line 12 performs the elitism process, which replaces E r of the chromosomes in the new population with the same fraction of best performing chromosomes in the old population. In line 15, if the fitness value of the best chromosome remains unchanged for a given number of consecutive iterations Q, the algorithm breaks out of the for loop and returns the best/optimal user association solution µ * k . Line 17 returns the optimal user association, which is the chromosome in the population with the best fitness value. This solution indicates which ANs the users should be associated with but does not give the number of BBUs that should be allocated to the users to meet their QoS requirements. Therefore, line 18 of Algorithm 1 inputs the optimal user association decision into Algorithm 2, and this returns a list Assoc i o j containing the AN o j (which is the gene µ k o j ,i ) serving user i, and the number of BBUs N BBU i determined according to (21) that are allocated to the user.

C. THE OPTIMAL SOLUTION
The optimal solution to the problem in (18) is obtained to validate the optimality of the proposed GA using the Gurobi solver. However, it is important to note the non-linearity of problem (18) introduced by the product of the variables µ o j ,i and ω o j ,i,c j in the objective function and constraints. Such a non-linearity hinders our usage of the Gurobi solver. To resolve this issue, the problem is linearised, as detailed in Appendix A. The reformulated problem in (28) is an integer linear programming (ILP) problem whose optimal solution can now be derived using the Gurobi solver via linear programming (LP) relaxation, Branch and Bound (BnB), and other advanced mixed integer programming techniques [63]. Select two parents from M using the roulette wheel technique 8: Carry out crossover on the two parents to create two children Break the for loop if fitness value of best chromosome does not change for Q consecutive iterations 16: end for 17: return the user association solution µ * Consequently, we utilize the Gurobi solver to obtain the optimal solution to the problem defined in (28).

D. COMPLEXITY ANALYSIS
In this section, we analyze the time complexity of the proposed GA versus the optimal solution, the baseline association, and the random user association used as benchmark solutions. The big Omicron (big-O) is employed to characterize the time complexity of the algorithms. The big-O is a mathematical notation that gives a measure of an algorithm's worst-case execution time or required memory in relation to the problem size. A detailed description of the big-O can be found in [64] and [65].

1) COMPUTATIONAL COMPLEXITY OF THE PROPOSED GA
The GA performs the selection, crossover, and mutation operators in each generation. Similar to many roulette wheel

2) COMPUTATIONAL COMPLEXITY OF THE OPTIMAL SOLUTION
The problem defined by (28) represents a variant of the well-known knapsack problem, which is NP-complete [67]. The optimal solution to such a problem can be obtained via LP relaxation and BnB, but this requires exponential upper bound time complexity in tandem with the problem size [67], [68]. Since the decision variables are binary, the problem's search space has a size of two to the power of the number of binary variables. Therefore, the time complexity of the optimal solution is given by O(2 |I|×(|B|×|C B |+|U |×|C U |+|H|×|C H |+|S|×|C S |) ), where |B|, |U|, |H|, |S|, denote the number of ANs in the MBS, LAP, HAP, and SatComs RAN respectively. On the other hand, |C B |, |C U |, |C H |, and, |C S | represent the number of BBUs owned by an AN in the MBS, LAP, HAP, and SAT-COMs RAN respectively. Since the worst-case time complexity of the optimal solution is exponential, algorithms that yield near-optimal solutions but with polynomial complexity should be considered; hence the GA proposed in this work.

3) COMPUTATIONAL COMPLEXITY OF THE BASELINE AND RANDOM USER ASSOCIATION
In this work, we also analyze the baseline and random user association (RUA) schemes as benchmark solutions.
The baseline association, also referred to as the greedy algorithm in this work, associates users with ANs based on maximum SINR. The description of the greedy algorithm, together with its pseudo code, is given in our work [32]. The time complexity for the computation of the SINR values from users to all ANs is of the order O(|I| × |O|) where |O| is the total number of ANs in the ITNTN. For each user, the SINR values to the different ANs must be sorted such that if the AN with the highest SINR does not have sufficient resources to meet the user's QoS requirements, then the user is associated with the next best AN. The time complexity due to sorting is O(|O| log |O| × |I|). Therefore, the overall complexity of the greedy algorithm is given as O(|I| × |O| + |O| log |O| × |I|) which can be reduced to O(|I| × |O| × log |O|). As the population size M and the number of generations G of the GA are usually greater than the number of access nodes |O|, it is clear that the greedy algorithm has a much shorter worst-case running time than the proposed GA. However, as the results will show in section V, the GA achieves a better performance as far as maximizing the objective function in (18) is concerned.
On the other hand, the RUA approach associates users randomly with any available AN. Such an algorithm has the shortest worst-case running time, which is proportional to the number of users. Hence, the time complexity of the RUA algorithm is given by O(|I|). Important to note is that, like the GA, both the greedy and RUA algorithms prioritize (1) the euRLLC use-case and (2) the use of NTNs over the TNs for service provisioning of the mobile LDHMC users.

V. RESULTS AND PERFORMANCE EVALUATION
In this section, the performance of the proposed GA is compared to the optimal, the greedy [32] and random user association (RUA) solutions. First, the main simulation parameters are presented, and after, the results and their discussion.

A. SIMULATION ASSUMPTIONS
We consider a circular urban region of 3 km radius that is within the coverage of a LEO satellite and a HAP AN and  also contains 1 LAP AN and MBSs with a radius of 2 km and 1 km respectively. Fig (5) is an example of the network deployment with 2 MBSs in the considered user distribution area. Users within the considered region are uniformly and randomly distributed. Table 2 [9], [24], [43], [44] gives the radio environment parameters used to validate the proposed solution. Given that the performance of the GA is highly dependent on the probability of crossover P c , probability of mutation P m and population size M, the values of these parameters are determined before results analysis.

B. GA PARAMETER SETTING
The appropriate parameters used in the proposed GA solution are identified in this sub-section. Fig. 6 shows the effect of P c on the GA convergence. P m , M and the number of users |I| are fixed at 0.1, 50 and, 80 respectively, while P c is varied from 0.2 to 1, with increments of 0.2. It is observed that the higher the value of P c , the larger the fitness value; hence, the better the solution found by the GA is at satisfying the objective function. Since P c = 0.8 performed well as per Fig. 6, it is chosen to be used in this work. Next, the effect of P m on convergence is analysed by setting P c = 0.8, M = 50 and |I| = 80. Fig. 7 shows that the higher the value of P m , the worse the performance of the GA, as the algorithm is transformed into random search. P m is set to 0.1 since its fitness value and convergence rate is much better than any other value of P m as can be observed in Fig. 7. In Fig. 8, the effect of population size M on the GA convergence is analysed, with P c = 0.8, P m = 0.1 and |I| = 80. The figure shows that convergence speed increases with the population size M . Also, we observe that convergence is achieved by the 200 th iteration for all population sizes. In the following section, we set the population size M to 50 and the number of iterations G to 150. These parameters are chosen to strike a balance between the accuracy and computational complexity of the GA, as both increase with M and G. Table 3 gives the parameters used for the GA.

C. SIMULATION RESULTS
To validate the performance of the proposed GA, we simulate the optimal solution based on the gurobi solver of the problem VOLUME 10, 2022  described in (28). In addition, we compare the proposed GA with the greedy algorithm [32] and the random user association (RUA) scheme.
Performance evaluation is based on three main metrics; the acceptance ratio (AR), the spectrum efficiency (SE), and the handoff probability. In this work, the user AR quantifies the ratio of served users to the total number of users in the network. On the other hand, the SE is the ratio of the overall system data rate to the total network bandwidth [69], while we define the probability of handoff as the ratio of the number of users that experienced a handoff (and are thus served by another AN) to the total number of mobile users during a given TTI.

1) IMPACT OF TRADE-OFF FACTOR α
First, we analyze the effect of the trade-off term α on data rate maximization (objective 1) and handoff minimization (objective 2) in (18). As α affects only mobile users, for this analysis, we consider 20 LDHMC users in a network comprising of 5 ANs, i.e., 2 MBSs, 1 LAP, 1 HAP, and 1 SatComs AN. Fig. 9 depicts that in solving the MOOP in (17) as a SOOP in (18) for varying values of α, a set of Pareto-optimal solutions exist. These solutions are generated using Algorithm 1 for values of α ranging from 0 to 1 with an increment size of approximately 0.0526. As Fig 9 shows, the generated Pareto-optimal solutions form a Pareto-optimal front below which the region comprises suboptimal solutions, and above which are infeasible solutions. From the figure, the points are concentrated at both ends of the curve. This shows that for the mobile users, α acts in a manner as to either maximize data  rate or minimize mobility-induced handoff. When the value of α is lower than 0.5, function two of (18) is maximized, which ultimately minimizes the handoff probability, and as α increases beyond 0.5, then function one is maximized consequently maximizing data rate. This observation is further supported by Figs. 10 and 11. In Fig. 10, the data rate is low for low values of α, and a step to higher values of data rate is observed at α ≈ 0.5. In the same manner, in Fig. 11, the probability of handoff is approximately 0 for α < 0.5 when objective two is prioritized, and once α increases above 0.5, the handoff probability increases since the priority becomes data rate maximization, and the nodes that maximize data rate are not necessarily the same as those that minimize mobility-induced handoff. The instability observed in both Figs. 10 and 11 for α > 0.5 is caused by motion of the users. The users keep moving at different velocities out and into coverage of different ANs resulting in unstable achieved total data rate and handoffs experienced. Since α either maximizes the data rate or minimizes the mobility-induced handoff, in all the following simulations, we assume that the objective of the mobile users is to minimize the handoff probability and, as such, set α = 0 for the LDHMC service group.

2) IMPACT OF USER DENSITY
We then evaluate the proposed algorithm's performance while varying the number of users in the network. We maintain the number of ANs at 5, with 2 MBSs, 1 LAP, 1 HAP, and 1 Sat-Coms AN.
In Fig. 12, we analyze the user AR performance of all algorithms. Generally, as the number of users in the network increases, the AR reduces due to resource scarcity. On average, the AR achieved by the GA, optimal, greedy, and RUA algorithms is 0.87, 0.86, 0.85, and 0.84, respectively. The proposed GA achieves an AR that is better than the optimal, greedy, and RUA solutions by 0.71%, 2.02%, and 2.75%, respectively. The GA performs better because, unlike the greedy and RUA algorithms, it is optimized to consider all different user association possibilities, thereby serving users with fewer ANs within their coverage first. Also, since the fitness value of the GA increases with the number of users admitted to the network, the proposed algorithm performs slightly better than the optimal algorithm that focuses on maximizing the data rate without regard to the AR. The greedy and RUA algorithms have more or less the same performance since we consider a network with a small number of nodes. Therefore, there is a high chance of selecting the same node within a user's coverage, whether by random or through the use of maximum SINR. Fig. 13 depicts the performance of the different algorithms with respect to SE. As the number of users in the network increases, the total data rate increases, thus increasing the achieved SE. However, at about 60 users and above, the SE remains constant since the available resources become insufficient to meet the QoS requirements of all users. Ultimately, all algorithms saturate, as the achieved total network data rate remains constant irrespective of the number of users in the network. The SE achieved on average by the GA, optimal, greedy, and RUA is 10.79, 10.79, 10.28, and 10.25 b/s/Hz respectively. The performance of the GA closely follows that of the optimal solution, outperforming the greedy and RUA algorithms by 4.81% and 5.094% on average, respectively. It is important to note that the RANs in the ITNTN have BBUs of different bandwidths. Therefore, a RAN may have a  higher SINR, but because of a smaller sized BBU, it achieves less data rate than another RAN with a bigger sized BBU. Therefore, an association based on maximum SINR does not guarantee the maximum data rate in the ITNTN, thereby leading to a lower SE when compared to the proposed GA, whose value function is based on maximizing the achieved data rate. The performance of the RUA is also lower than the GA, as this algorithm associates users randomly to any available capable RAN without any regard for the achieved data rate. The excellent performance of both the GA and the optimal solutions is because both these solutions are based on the maximization of the data rate of the ITNTN.
In Fig. 14, the euRLLC user acceptance ratio performance is depicted. The euRLLC users are prioritized over other users for all four algorithms, as it is vital to mitigate call blocks for this use-case. Consequently, for the number of users from 10 to 70, all euRLLC users are accepted into the network since resources are still available to meet their needs. For the number of users beyond 70, a fraction of euRLLC users is dropped due to limited resources to serve all missioncritical users. The GA and optimal algorithms have the same performance with an average AR of 0.99, performing better than the greedy and RUA, with average euRLLC AR of 0.98 and 0.97, respectively. The GA and optimal solutions perform better than the greedy and RUA due to their intelligence in considering all the options necessary to mitigate euRLLC user call blocks. The large coverage NTNs were prioritized in all four algorithms to serve the mobile LDHMC service group. Therefore, as is depicted by Fig. 15, the AR for this service group is 1 for all algorithms and all number of users. This is because of available resources in the network to serve this service class. On the other hand, Fig. 16 shows the acceptance ratio of the feMBB use-case. In this work, the priority of this use-case is lower than other use-cases, hence the steep decline in AR beyond 50 users when the resources in the network are no longer enough to serve all users. Because of its intelligence in considering the different available AN association options, the GA has a better feMBB of 0.79 on average than the optimal, greedy, and RUA, which have an average of 0.78, 0.76, and 0.76, respectively. We analyze the handoff performance of all algorithms in Fig. 17. The optimal and GA algorithms can associate the LDHMC mobile users to the ANs with the largest cell radius and thus achieve a handoff probability of 0 for all numbers of users in the network. On the other hand, the handoff probability for the greedy and RUA algorithms is worse by 16.7% and  41.3% on average, respectively. This is because the LDHMC user association keeps alternating between the NTN ANs depending on the maximum SINR for the greedy algorithm and randomly for the RAU algorithm. The handoff probability reduces with the number of users for these two algorithms because resources of smaller radius ANs are depleted, and users are now forced to associate with the large coverage cells.

3) IMPACT OF ACCESS NODES DENSITY
Next, we analyze the impact of AN density on the four different algorithms. We vary the number of MBS from 1 to 6 while maintaining the number of LAPs, HAPs, and Sat-Coms to 1 AN each. The number of users in the network is maintained at 80. Fig. 18 shows that as the number of ANs in the network increases, the SE also increases since resources to support the users' data rate requirements keep increasing. The SE performance of the GA is within an average of 0.4% of the optimal SE and outperforms the greedy and RUA algorithms by 1.23% and 0.97% on average, respectively. Fig. 19 shows that the GA outperforms all the other three algorithms in terms of AR, with an average AR of 0.867 compared to 0.855, 0.774, 0.801 of the optimal, greedy, and RUA, respectively, translating to 1.41%, 10.8%, and 7.6% better performance respectively. Since all users in the network cooperatively contribute to the fitness value of the GA algorithm, this solution maximizes the number of admitted users in the network and hence achieves a higher AR. The greedy algorithm has the worst performance in terms of SE and AR since the association is based on maximum SINR without regard for the user and AN distribution, and hence lacks the intelligence of first associating users that are within the coverage of few ANs.  Moreover, Fig. 20 also shows that as the number of ANs in the network increases, the GA achieves superior performance in terms of handoff probability. The handoff performance of the optimal algorithm falls below that of the GA by 8.4% on average since this algorithm chooses to maximize the total data rate in the network at the expense of minimizing handoff. On the other hand, the GA performs better than the greedy and RUA in terms of handoff probability by 14.9% and 51.8% on average, respectively. This is because the greedy chooses an NTN AN with the highest SINR in each TTI, while the RUA chooses any NTN AN at random. Consequently, the ANs chosen for the association by the greedy and RUA algorithms   keep changing in each iteration, while the GA fitness value is formulated so that the mobile LDHMC users are associated with the AN with the largest cell radius in each TTI. These results demonstrate that the proposed GA is well suited for future scenarios characterized by highly mobile users for which increased handoff implies increased delays and a high probability of call drops due to handoff failure, consequently degrading the QoS of the users.

4) IMPACT OF NETWORK OVERLOAD
In the previous simulations, the argument is that the nodes with the broadest coverage should be prioritized to serve the LDHMC service group to reduce mobility-induced handoff. However, in this section, we analyze the performance of the proposed algorithm in a network experiencing overload  conditions. To best evaluate this scenario, we decided to give the same priority to both the feMBB and LDHMC service groups so as to have a fair comparison in overloading conditions for both use-cases. The euRLLC use-case is still prioritized, and the simulation is performed considering 5 ANs, i.e., 2 MBS, 1 LAP, 1 HAP, and 1 Satcoms AN. The objective is to analyze how the system responds to the distribution of resources to the different use-cases.
First, Fig. 21 shows that for all algorithms, there is a small change in the SE as the number of users increases in the overloading conditions. In this condition, the limit on the number of users the network can support has already been reached, beyond which the network observes only a slight variation in SE. Nonetheless, on average, the GA still outperforms the greedy and RUA algorithms by 0.97% and 0.9% and is within 0.5% of the optimal SE. The Figs. 22, 23, and 24 depict the AR of the euRllC, LDHMC, and feMBB use-cases respectively. In all figures, the AR of the respective use-cases decreases as the number of users keeps increasing beyond the threshold that the network can support. Fig. 25 is a combination of Figs. 22, 23, and 24, showing the average AR of the different algorithms, for the different use-cases.
As can be observed in Fig. 25, All algorithms prioritize the euRLLC use-case over the other two, with the GA being better than the optimal, greedy, and random algorithms by 3.5%, 8.55%, and 5.68% respectively. This continues to show the strength of the formulated fitness function of the GA in prioritizing mission-critical users.
Since the LDHMC and feMBB have the same priority, we observe an almost similar AR for the two use-cases for all algorithms in Fig. 25. Fig 26 shows the AR of all algorithms in overloading conditions. The GA is able to still outperform the optimal, greedy, and RUA algorithms by 0.37%, 7%, and 4.78%, respectively, and yet still maintain a low mobility-induced handoff probability as observed by Fig. 27.

VI. CONCLUSION
In this paper, we have formulated a user association and resource allocation problem in the ITNTN as a MOOP that maximized the total network data rate and minimized the mobility-induced handoff. Moreover, the mission-critical euRLLC service group provisioning was prioritized over other service groups. The MOOP was transformed into a weighted sum SOOP, which was solved using the GA. In the GA, service group-dependent fitness functions were formulated to determine the near-optimal user association and resource allocation solution. The Simulation results showed that for an increasing number of ANs, the proposed GA's SE is within 0.4% of the optimal solution. At the same time, the GA's user AR and handoff probability outperformed the optimal, the greedy SINR association-based, and the random user association solutions. While the greedy and RUA algorithms are characterized by a shorter running time compared to the GA, the above results show that the GA achieves a better SE and lower probability of handoff. In future work, we plan to investigate the proposed service-aware user association and resource allocation in the ITNTN based on reinforcement learning (RL). In RL, once the training is done, the agent can make decisions in real-time, a considerable advantage for wireless networks, especially those serving mission-critical users.

APPENDIX A
The multiplication of the two decision variables µ o j ,i and ω o j ,i,c j introduces a non-linearity in the objective function and constraints of the optimization problem formulated in (18). Similar to [70], a linearisation term is introduced that replaces the product with a single binary variable as depicted by (24) to avoid this non-linearity. o j ,i,c j = µ o j ,i ω o j ,i,c j , ∀ o j ∈ O, ∀ i ∈ I, ∀ c j ∈ C j (24) o j ,i,c j ∈ {0, 1} in (24) is a binary decision variable that is one when a user i is associated with AN o j and allocated a BBU c j , while it is zero otherwise. Subsequently, the product µ o j ,i ω o j ,i,c j in the objective function and constraints of (18) is replaced with o j ,i,c j . Furthermore, the linearised optimization problem must include additional constraints that establish the relationship between o j ,i,c j , µ o j ,i , and ω o j ,i,c j ; defined as o j ,i,c j µ o j ,i , ∀ o j ∈ O, ∀ i ∈ I, ∀ c j ∈ C j (25) o j ,i,c j Henceforth, the modified integer linear programming (ILP) problem is formulated as s.t. Important to note is that the data rate L o j ,i,c j achieved via a BBU c j by a user i associated to AN o j is a constant since this term is computed before hand. Also the terms π o j ,i and ρ iυ are known before hand. He is currently a Senior Lecturer with the Department of Electrical Engineering, University of Cape Town. His research interests include system security, machine learning, and wireless sensor networks, primarily focusing on the IoT applications and cyber-physical systems.
NECO VENTURA (Life Member, IEEE) is currently a Senior Research Scholar, the Head of the Centre for Broadband Networks, and the Director of the Communications Research Group, Department of Electrical Engineering and the Department of Computer Science, University of Cape Town (UCT). He has held positions on several conference organizing committees. He is on the technical program committees of various international conferences. Over the last decade, he has contributed over 100 publications in refereed journals, chapters in books, and refereed international conferences. His research interests include the field of networking, currently centered on next generation mobile networks and architectures, the Internet of Things, machine to machine communications, SDN, NFV fog/edge computing, and 5G. He is a member of the IEEE Computer and the IEEE Communications Societies. VOLUME 10, 2022