Design and Optimization of a UAV-Enabled Non-Orthogonal Multiple Access Edge Computing IoT System

Unmanned aerial vehicle (UAV)-enabled mobile edge computing (MEC) in IoT systems has recently gained prominence as a solution to accommodate the low latency and high energy efficiency needs of emerging and developing 5G-and-beyond IoT applications. Given the limited spectrum available, non-orthogonal multiple access (NOMA) communication scheme provides a spectral- and power-efficient task offloading means. In this work, a novel framework is proposed for a multiple-UAV-enabled MEC IoT system incorporating uplink NOMA. The goal is to maximize the sum bits offloaded from all the IoT devices (IoTDs) in a finite service period. This is achieved through joint optimization of the associations between the IoTDs and the UAVs, the IoTD transmit powers, as well as the UAV trajectories, given UAV and IoTD energy budgets and quality-of-service (QoS) criteria. Towards that end, a specialized penalty block coordinate descent (P-BCD) algorithm is proposed to decompose the original problem into 4 subproblems that are solved alternately and iteratively. By direct comparison with several benchmark schemes including orthogonal multiple access (OMA), our proposed scheme displays performance enhancement in terms of sum bits offloaded and energy expenditure. Finally, a performance comparison is carried out for the proposed scheme with different number of employed UAVs to find the most cost-effective number of UAVs to deploy.


I. INTRODUCTION
The past decade has been a period of proliferation and evolution of the Internet of Things (IoT) technologies. The estimate that by the year 2025, more than 30 billion IoT devices will be in use is becoming within reach [1]. Applications pertaining to the fields of augmented and virtual reality (AR/VR), autonomous driving, smart cities, agriculture, disaster management, search and rescue operations and real-time video analytics are in high use and demand [2], [3], [4], [5]. Such applications are typically enabled by the exchange of high data traffic between a myriad of computationally-capable and environmentally-aware networked devices. A simple The associate editor coordinating the review of this manuscript and approving it for publication was Xiaolong Li . example would be the integration of multiple sensors, meters, and street and traffic lights that communicate with each other to enhance city infrastructure and reduce energy consumption and carbon footprints.
In addition to the above, the recent miniaturization of wireless transceivers has made it possible to carry a complete base station (BS) onto an unmanned aerial vehicle (UAV), making them not only capable of wireless communication as users, but also as Aerial Base Stations (ABSs). A clear advantage of employing UAV-ABSs to serve ground terminals is the ability to readily establish line-of-sight (LoS) air-to-ground (A2G) links due to the high altitudes of UAVs [6]; this then mitigates potential signal blockage and shadowing, and hence enhances the channel strengths for better capacities. The low-cost of UAVs also means that dynamic changes in networks and 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/ service requests can be more effectively handled through carefully adapting UAVs deployment. Consequently, UAV-ABSs have found widespread use in many 5G and beyond-5G (B5G) IoT networks, especially massive IoT scenarios. This is due to their ability to hover close to IoT devices to establish higher-rate LoS links, effectively lowering the energy consumption, and to relay the offloaded task input data to the nearest ground base station to be processed by servers at the centralized cloud of the network. However, to accommodate the increasingly time-sensitive requirements of recent and emerging IoT applications, data and services frontiers can be pushed away from the centralized cloud to the edge of the network, allowing task computation and processing to be held in proximity to the IoT devices thereby reducing latency while maintaining QoS. This is typically referred to as mobile edge computing (MEC) [7]. Specifically, in many recent IoT scenarios, MEC provision has been implemented by mounting cloudlets (or edge servers) on the UAVs for computational capabilities, promising better servicing times [8]. UAV-assisted MEC comes with its own set of challenges that need to be overcome in order to tap into the benefits. First, UAVs have limited onboard energy, and to accommodate that fact, one should aim to optimize the energy consumed through on-board computation in the designated service time. This should jointly be done while planning energy-aware trajectories for the UAVs to reduce wasted energy due to flight as a second challenge. Third, communication resource allocation at the UAVs should be made to support all the latency-sensitive services in the IoT network, while also keeping in mind the smaller contribution (compared to UAV flight) to the energy efficiency of the system [8], [9]. It is important to note that the three challenges mentioned earlier are all coupled such that proper design of an energy-efficient UAV-MEC IoT system to satisfy all the QoS and service time requirements of the IoT devices while considering the energy budget of the UAVs is a highly complex problem [9]. Furthermore, while simpler to plan and accommodate, a single UAV is usually not capable enough to handle all the computational load from the devices and may not have the time to hover over all the users for as much as needed to offload and compute all the task input data within the time restriction. Hence, a practical solution is to deploy multiple UAVs that jointly operate to satisfy the system requirements [10].
Another challenge that MEC commonly faces is that of limited spectrum for offloading task data (and to a lesser extent retrieving the computation output due to its significantly smaller size). This could mean that traditional communication schemes that allocate time and/or frequency resources among devices may fall short of supporting higherend IoT applications. The simple option of increasing transmit powers of the devices to increase the offloading rate is not wise, as this could result in interference with other neighboring networks as well as cutting down on the IoT devices longevity. A possible solution is to employ nonorthogonal multiple access (NOMA) [11] in the uplink, thus enabling more efficient task offloading management.
In uplink NOMA, superposition coding (SC) is first employed to perform power-domain multiplexing at the devices end. Then, on the UAVs end, successive interference cancellation (SIC) is performed to demultiplex the desired signals [12]. Typically, the signals are decoded in descending order of channel gains of the device-UAV links. Through NOMA, multiple users can share the same limited time and/or frequency resources with improved spectral efficiency.
Motivated by the above discussions, in this paper, we study and tackle an energy-efficient multiple-UAV-assisted MEC IoT framework incorporating uplink NOMA, which by itself is a novel framework that needs more exploration. To the best of our knowledge, this is the first paper of its kind to jointly do the following: 1) Enable the UAVs to dynamically adjust their locations, incorporating the resultant trajectories as part of the optimization procedure. (as opposed to the very common device clustering for NOMA using heuristics and confining the UAVs to the designated cluster centers); 2) Utilize the penalty method in conjunction with the block coordinate descent (BCD) method for optimizing the UAV-device dynamic binary associations. (as opposed to the common relaxation of the binary association variable constraints in such cases, which reduces the accuracy); 3) Perform a direct comparison between the NOMA and orthogonal multiple access (OMA) counterparts of the proposed framework; 4) Study the effect of number of UAVs employed on the system performance.
The rest of the paper is organized as follows: Section II provides a literature review of multiple-UAV-assisted MEC networks followed by outlining the major contributions of this paper. Section III presents the system model and attributes. Section IV deals with the problem formulation and proposed solution method. Section V discusses the simulation results and numerical analysis. Finally, Section VI concludes the paper and presents possible future work.

II. RELATED WORKS AND NOVELTY OF THE PROPOSED WORK
In this section, we present a summary of the related works in the literature that studied UAV-assisted IoT systems and proposed solutions to various related optimization problems and highlight our contributions.
In the first category of relevant research works on UAV-assisted IoT systems, the focus simply lies on the collection of sufficient (possibly time-sensitive) data from IoT devices to meet certain QoS or latency demands, with the computing process neglected. For example, in [13], Mozaffari et al. studied a multiple-UAV-assisted IoT framework dealing with a time-variant device activation process, minimizing the total transmit power of the IoT devices by jointly optimizing the dynamic UAV locations and flight paths, the device-UAV associations, and the device powers subject to devicespecific signal-to-interference ratio (SINR) constraints. Also, Samir et al. [14] investigated the practical problem of deadlines on data collection from IoT devices. The authors aimed to optimize the trajectory of a UAV ABS in collecting data from the maximum number of served devices, where a minimum amount of data is to be uploaded from a device for it to be considered served. Furthermore, Yang et al. [15] investigated a multi-UAV-enabled data dissemination system serving a set of IoT nodes with unsaturated data requirements spread across a post-disaster region. The authors aimed to optimize the mission completion time and UAV transmit powers and trajectories. They utilized the BCD method in combination with the bisection search method to obtain a suboptimal solution to the original intractable problem. In an effort to overcome the inherent access delay in multi-user UAV-enabled employing OMA, Lu et al. [16] proposed a UAV-enabled system incorporating uplink NOMA for data collection. The authors aimed to maximize the users' sum rate by jointly optimizing the deployed UAV position and the transmit powers of the users. The authors show notable enhancement of the sum-rate over conventional OMA and terrestrial NOMA. Duan et al. [17] studied a multiple-UAV-assisted IoT system incorporating NOMA scheme in the uplink. The authors aimed to maximize the sum of uplink achievable rates by jointly optimizing the devicesubchannel assignment, the device transmit powers, as well as the UAV altitudes. Better convergence behavior as well as sum-achievable-rates were obtained compared with the OMA scheme.
In the second set of works, edge computing (in particular, MEC) plays an important role in the UAV-assisted IoT systems, shifting the focus to the more useful computed/processed data size as a performance metric, in addition to making good use of the computational resources of the edge computing servers. In [18], Zhang et al. investigated a UAV-assisted MEC IoT system in which a single UAV can act as both a computation assistant to the IoT devices as well as a decode-and-forward (DF) relay for transferring offloaded bits from the devices to a more computationally capable access point (AP). The authors aimed to minimize the total energy consumed in the system through communication, computation, and UAV flight by optimizing the computational resource allocation, TDMA time slot scheduling, power allocation, and the UAV trajectory. Along similar lines, Yu et al. [7] studied a UAV-assisted MEC IoT system where a single UAV is used to provide computational services as well as facilitate the exchange of information between the devices and several terrestrial edge clouds (ECs), with the aim of minimizing the weighted sum of the service delay of the IoT devices and the total UAV energy consumption. Also, Zhang et al. [19] studied a multiple-UAV-MEC IoT system employing the partial offloading computation mode, making use of the IoT devices innate computational capability to aid the process. The authors aimed to maximize the number of served IoT devices through joint optimization of the service indicator variables, the time and computation resource allocation scheme, as well as the UAV trajectories. They utilized the BCD method and successive convex approximation (SCA) technique for obtaining a sub-optimal solution. In an effort to meet more sensitive access delay needs of users, Diao et al. [20] studied a UAV-enabled uplink-NOMA-based MEC system. The authors focused on enhancing energy savings and user fairness through minimizing the largest energy consumption among users. This is done by joint optimization of trajectory, as well as communication and computation resource allocations. The authors showed considerable energy savings as well as lower energy-consumption gap across the served users compared to benchmarks. Due to its efficiency, NOMA was incorporated into the multiple-UAV-enabled MEC IoT system studied by [21] where the authors aimed to maximize the total energy efficiency while minimizing the number of UAVs and the total cost of UAV and central cloud use. This is done through jointly optimizing the transmit powers, computation and bandwidth allocation, UAVs 3D locations, and designation of the UAVs to the NOMA device clusters. They utilized the BCD method for handling the problem. Similarly, Zhang et al. [22] were concerned with establishing a multiple-UAV-MEC IoT network incorporating uplink NOMA with the goal of minimizing the overall weighted cost of energy consumption in the system by joint optimization of radio and computation resources allocation and UAV trajectories. The authors demonstrated the effectiveness of the proposed scheme in meeting its energy efficiency goals, outperforming the OMA counterpart with greater improvements as the size of the system grows. Tackling the issue of secure UAV communications provision, Xu et al. [23] proposed UAV-assisted MEC frameworks for TDMA and NOMA schemes. In a setting with multiple available eavesdroppers, the authors utilize two UAVs, one for service provision, and the other for jamming and hence protection against eavesdropping. The authors aim to maximize secure computing capacities by jointly optimizing computation resources, communications resources, and UAVs' trajectories. They adopt a P-BCD algorithm to obtain a suboptimal solution to the original problem. The authors show that the NOMA scheme outperforms its TDMA counterpart in terms of security capacity. A similar model is proposed by Lu et al. [24]; a NOMA-based dual-UAV-MEC system with a flying service-provision UAV, a flying eavesdropping UAV, devices to be served, and a terrestrial jamming device. The authors formulate an optimization problem maximizing the average system security computation capacity by jointly optimizing SIC ordering variables, transmit powers, locally computed bits, UAV-MEC CPU frequency, and the service-providing UAV trajectory. The maximum error in estimating the motion of the flying eavesdropper is taken into the account. The authors illustrate the need for joint optimization of all the included variables to achieve the best secure computation capacity compared with several benchmarks. Despite the popular approach of energy consumption minimization in many of the reviewed works on UAV-assisted MEC systems for yielding notably higher efficiency in comparison to traditional schemes and benchmarks, in the system model that is proposed next, we take a different approach by assuming that a means of recharging or replenishing the devices is more accessible, reducing the importance of system energy consumption minimization. In contrast, we assume that the IoTDs lack the ability to compute a portion of the task-input data locally, relying exclusively on offloading their tasks for getting meaningful output. Hence, to ensure that we compute as much of the task-input data as possible within the limited time frame, we focus on extracting as much data from the IoTDs as possible for given operating period and energy budgets.
Based on the discussions above, the contributions of this work can thus be summarized as follows: • We propose a novel framework of a multiple-UAVenabled MEC IoT system utilizing NOMA for uplink task-input data transmissions, where edge computing services are provided to a collection of IoTDs. The UAVs cooperatively serve the IoTDs to meet their QoS demands for offloading data to be computed at the UAVs. The UAVs are to dynamically adjust their locations to increase the flexibility and efficiency of edge computing service provision. In our system, the available bandwidth is partitioned equally among the UAVs; hence, no co-channel interference is present.
• We formulate an optimization problem with the objective of maximizing the sum bits offloaded from all IoTDs by jointly optimizing the UAV-IoTD associations, the IoTD transmit powers, as well as the UAV trajectories. Our formulation takes into account the energy budgets of both the IoTDs and the UAVs, as well as the QoS requirements representing the minimum quantity of bits that need to be offloaded from the IoTDs for useful operation. The formulated optimization problem is a mixed-integer non-linear programming (MINLP) that is challenging to solve.
• We propose a penalized block coordinate descent (P-BCD) based algorithm to deal with the binary association optimization variables, where four steps are iteratively executed for variable-block optimization. The SCA technique is applied in the second to fourth steps to tackle the non-convex objective and constraints. All UAVs considered are assumed to carry an onboard communication circuit in addition to their onboard computing processor, enabling data exchange and task computation to occur simultaneously. All devices in the system are assumed to carry a single antenna. Let T denote the task-offloading period, during which the UAV-MECs fly from a starting location and return back at the end of the period. For ease of exposition, the computing cycle period T is discretized into N equal-length time slots of length δ t , where N denotes the set of time slots. The discretization is made such that the slot length δ t is small enough to consider the UAVs as being static within each time slot.

B. COORDINATE SYSTEM AND COMMUNICATION CHANNEL MODELING
We consider a 3D Cartesian coordinate system for spatial analysis where the fixed horizontal location of IoTD k ∈ K is denoted by the vector w k ∈ R 2×1 and all IoTD locations are assumed to be known to all UAVs beforehand. The location of UAV m ∈ M at time slot n ∈ N projected on the ground plane is denoted by the vector q m [n] ∈ R 2×1 and the UAVs are assumed to fly at a fixed altitude H where the periodic flight trajectory of UAV m is denoted by the ]. The average speed of UAV m within time slot n is denoted by v m [n] and is given by These speeds are all upper bounded by a maximum The UAV trajectories need to satisfy the following constraint: implying that each UAV is to return to its starting location after the period of N time slots has passed to serve the IoTDs in a cyclic manner. Further, another set of constraints is imposed to prevent collision between UAVs, Here, d min denotes the minimum safety distance. Given the open rural setting of the multiple-UAV-enabled IoT system, the influence of terrains, obstructions and shadowing effects can be ignored. This is confirmed by recent field experiments, which have proven that, even in an urban environment, the channels between the UAVs and the IoTDs are approximately modelled as line of sight [6]. Furthermore, the Doppler effect resulting from UAV mobility is assumed to be well compensated for at the receiver side [9], [25]. Hence, the channel power gain from IoTD k to UAV m in time slot n follows the free-space path loss model, and is expressed as where ρ 0 denotes the channel power gain at the reference distance of 1 m.

C. MULTIPLE ACCESS AND ACHIEVABLE RATES
We denote the total bandwidth pre-allocated to uplink data offloading from the IoTDs to the UAVs as B, which is partitioned equally between the UAVs, leading to a bandwidth of B/M per UAV. A binary association variable α m,k [n] is defined to indicate whether IoTD k offloads data to UAV m in time slot n, where a value of 1 indicates service and 0 no service. We adopt the common assumption that each IoTD can only offload data to at most one UAV in each time slot [10], which can be imposed mathematically as As mentioned earlier, power-domain NOMA is incorporated into the network and SIC is carried out on the superposed multiuser signal transmitted to an associated UAV. To this end, a ranking of the IoTDs in terms of channel gains is done for each UAV-MEC in each time slot to track changes in channel gains due to changes in UAV positions.
where N 0 denotes the power spectral density of additive white Gaussian noise at receivers.

D. ENERGY CONSUMPTION
In our model, downlink transfer of task-output data from the UAVs back to the IoTDs is neglected since the size of task-output data is much smaller than the size of task-input data and the download process hardly takes up much communication resources like bandwidth and power. Moreover, neglecting the download part reduces the model complexity and has been done by numerous studies, e.g., [22] and [26]. The total energy for task offloading is upper bounded by a known E k ; hence, To meet the QoS of IoTD k, the total number of bits offloaded in the period N , denoted by S k , is required to be greater than a preset minimum guarantee D k ; hence, Based on [27], the energy consumption for computing the offloaded task-input data from IoTD k at UAV m is modeled as VOLUME 10,2022 where κ U denotes the effective capacitance coefficient of IoTD k, which depends on its processor chip architecture [28], S m,k denotes the offloaded bits from IoTD k to UAV m in the period N , f CPU denotes the CPU frequency of UAV m during time slot n for computing the offloaded data of IoTD k, and C k denotes the number of CPU cycles required to process one bit of task of IoTD k. For simplicity, we approximate the UAV flight as a series of uniform rectilinear motions, each between two positions separated in time by δ t . This assumption is widely adopted in many related studies, including [26], [29], and [30]. Hence, no acceleration or turning radius is considered when modelling the UAV flight energy consumption. Let W m denote the mass of UAV m, then, based on the model detailed in [28], [31], and [32], the flight energy consumption for UAV m is given by Here, constraints (12b) reflect that the total UAV flight and computation energy consumption should not exceed the UAV energy reserves E 0 . Problem P 0 is non-convex due to the objective function (12a), the QoS constraints (9), and the UAV energy consumption constraints (12b), where the three sets of variables A, Q and P are highly coupled. Further, the binary nature of the IoTD-UAV association set A as well as the collision avoidance constraints in (4) impose two additional causes of non-convexity. Therefore, problem P 0 is an MINLP problem that is challenging to solve directly.

B. PROPOSED SOLUTION METHOD
To tackle the challenges proposed in the previous section for obtaining a suboptimal solution to P 0 , we propose a four-step P-BCD algorithm based on the penalized method [33], [34], [35], [36].  (13) is equivalent to the binary constraints (12c).
Next, we incorporate the equalities in (13) into the objective function to yield a penalized objective and hence, a penalized problem P 1 can be formulated as It can be proven that when λ is sufficiently large, problem P 1 is equivalent to problem P 0 [10]. Problem P 1 remains a nonconvex problem. To solve this problem, we consider applying the BCD method to optimize the decision variables in an inner loop and update the penalty factor in an outer loop until some convergence criterion is attained. Rather than optimizing the complete set of optimization variables in one problem, we optimize the four individual sets of variables A,Ã, Q and P individually keeping the other 3 sets of variables fixed. This is carried out iteratively until all the variables have been optimized. This is detailed as follows.

1) OPTIMIZINGÃ WITH FIXED A, Q, AND P
Upon fixing A, Q and P, problem P 1 becomes a minimization problem of the penalty term with respect toÃ with no constraints. Hence, the optimal value ofα m,k [n] at the rth iteration can be easily obtained by setting the first derivative of the objective with respect toα m,k [n] equal to zero. Hence, we get the following closed form solution: 2) OPTIMIZING A WITH FIXEDÃ, Q, and P Upon fixingÃ, Q and P at the rth iteration, problem P 1 can be formulated as Problem P 2 is non-convex due to the non-concave objective function and non-convex constraints (9) and (12b). We can solve this problem by applying the SCA technique [37], which can approximately solve non-convex problems by approximating them with convex problems in each iteration, followed by updating the values of the optimization variables iteratively. To tackle the non-concave S k term appearing in the objective and the LHS of (9), we rewrite the rate expression R m,k [n] as the difference of two concave log functions, i.e., R m, Hence, a concave lower bound approximation of S k in (9) is obtained such that and a convex upper bound approximation of S m,k is obtained: which can be plugged into (10) to get an upper-bound approximation E comp,α,ub,r m,k . Finally, problem P 2 can be approximated as In problem P 3 , the objective function is convex, and all constraints are convex. Hence, P 3 belongs to the class of standard convex optimization problems for which many standard convex optimization solution methods and solvers such as CVX [38] can be leveraged. We note that the feasible set for problem P 3 is always a subset of problem P 2 due to the lower bound approximations utilized. Thus, the optimal objective of problem P 3 serves as a lower bound of the optimal value of P 2 generally.
3) OPTIMIZING P WITH FIXEDÃ, A, and , Q Upon fixingÃ, A, and Q at the rth iteration, problem P 1 can be formulated as Problem P 4 is non-convex and the causes of non-convexity are the same as those encountered in the previous step, i.e., α m,k [n] and p k [n] always appearing together as a product in the sources of non-convexity. Hence, we follow the exact same procedure used before for deriving a concave lower-bound approximation forŜ asŜ p,lb,r and a convex upper-bound approximation for E comp m,k as E comp,p,ub,r m,k , where first-order Taylor expansion is applied with given local point P r = {p r k [n]} at the rth iteration on the same terms as before. The end result is that P 4 can be approximated as (24b) Problem P 5 is convex since the objective and all the constraints are convex; therefore, this problem can be practically solved using CVX or several other standard optimization techniques that are readily available.

4) OPTIMIZING Q WITH FIXEDÃ, A, and P
Upon fixingÃ, A, and P at the rth iteration, problem P 1 can be formulated as Problem P 6 is non-convex due to the non-concave objective function and the non-convex constraints (4), (9), and (12b). To tackle the non-concave S k term appearing in the objective and the LHS of constraints (9), we again rewrite the rate expression R m,k [n] as the difference of log functions, i.e., R m,  and problem P 6 is then reformulated as It can be shown that at the optimal solution of P 7 , constraints (25b) are all met with equality; otherwise, χ m,k [n] can always be increased to attain equality at no cost to the objective value assuming all other constraints and variables remain unchanged. Furthermore, constraints (25a) are obtained by squaring both sides of the respective constraints (4) for simpler derivation. Therefore, problem P 7 is equivalent to problem P 6 . P 7 remains a non-convex problem due to the non-concave objective function and non-convex constraints (25a)-(25d). Algorithm 1 Proposed BCD-Based Algorithm for Solving P 1 1: Initialization: Give feasible point for A,Ã, P, and Q.
Set iteration count r i = 0. Set a tolerance for solution accuracy i > 0. Specify the maximum number of iterations r i,max . 2: repeat 3: UpdateÃ via (14) in Step 1.

4:
Update A by solving P 3 in Step 2.

5:
Update P by solving P 5 in Step 3. Update r i = r i + 1 8: until the relative decrease of the objective is less than i or r i > r i,max We adopt the SCA technique to obtain a convex approximation of the problem that is then used to obtain the optimal values of the UAV trajectories by updating the solutions to the approximate problem iteratively. The convex terms q m [n] − w k 2 in the LHS of (25a) and (25b), respectively, can be tackled by applying the first-order Taylor expansion to derive a lower bound approximation of the two respective terms as follows, where the expansions hold at the rth iteration: Update the penalty factor λ r+1 = cλ r

5:
Update r = r + 1 6: until the relative decrease of the objective is less than or r > r max Next, to tackle the non-convexity of constraints (25c), we derive a lower-bound concave approximation of R m,k [n] by first-order Taylor expansion at the rth iteration, yieldinĝ Problem P 8 is a convex problem and many existing methods and solvers, such as CVX, can be used to solve it. For the same reasons explained in the previous steps, the feasible set for P 8 is a subset of that for P 7 and hence, its optimal objective serves as an upper bound to that for P 7 .

C. OVERALL P-BCD FRAMEWORK AND ANALYSIS
The proposed BCD-based algorithm for solving problem P 1 in the inner loop of the P-BCD framework is summarized in Algorithm 1. The overall P-BCD algorithm for solving P 0 is summarized in Algorithm 2. In the inner loops, problem P 1 with fixed penalty factor λ is approximately solved by applying the BCD-based Algorithm 1. In the outer loops, the penalty factor is updated to increase the penalty for violating the equality constraints. The convergence of the proposed P-BCD algorithm is verified and studied in-depth in [39], where it is proven using the Robinsons condition that the sequence of solutions generated from the P-BCD algorithm converge to a set of stationary solutions of the original problem. In practice, it is enough to rely on the relative change in the objective value of problem P 1 as an indicator of when to terminate the P-BCD algorithm.
The complexity of the proposed P-BCD algorithm (Algorithm 2) can be computed by first computing the complexity of the inner loop, which in turn depends on the complexities of the solutions to the subproblems in steps 3-7 of Algorithm 1. UpdatingÃ via (14) gives us a complexity of O(MKN ). For solving the standard convex non-linear optimization subproblems P3, P5, and P8; updating A, P, and Q, respectively; we utilize CVX with a primal-dual interiorpoint method of general complexity O(ς 3 log(1/ε)), where ς denotes the number of optimization variables and ε denotes the tolerable duality gap [41]. We assume that ε is constant for all 3 subproblems, and let I 1 and I 2 represent the number of iterations of the inner and outer loops, respectively. Hence, the overall computational complexity of Algorithm 2 can be found as O(I 2 I 1 (((MKN ) 3 + (KN ) 3 + (MN ) 3 ) log(1/ε) + MKN )). Hence, the proposed algorithm runs in polynomial time.

V. SIMULATION RESULTS AND DISCUSSIONS
In this section, we present numerical results to validate the effectiveness of our proposed scheme. We further investigate the impact of some system parameters on the performance of the overall system model.

A. TRAJECTORY INITIALIZATION SCHEME
In general, for iterative-based algorithms, such as our proposed Algorithm 2, both the converged solution and the system performance depend on the initialization scheme employed. The crux of the initialization procedure lies in how to obtain a good starting plan for the UAV trajectories with a general objective(s) that resembles those behaviors that can be inferred from our problem objective and formulation. We aim for initial trajectories that can essentially cover as much of the area containing the IoTDs as possible while being sufficiently separated and non-overlapping to achieve a more defined clustering of the IoTDs to get higher NOMA rates. We can get as close to achieving those objectives as possible by implementing simple circular trajectories in a circle packing scheme. To illustrate, we implement the initial trajectory for each UAV as a circular trajectory with a uniform VOLUME 10, 2022  . Next, the radius of the smallest circle with center c g that covers all IoTDs, denoted by r u , is determined, which is simply the maximum distance between the center c g and all IoTDs, i.e., r u = max k∈K w k − c g . We utilize the circle packing (CP) scheme detailed in [40] with the given number of UAVs employed M and enclosing circle radius r u to find the centers of the M packing circles, c m trj , and the common radius, r cp , of each circle. Rather than implementing trajectory circles of radius r cp , we halve the radius of each circle to balance the number of IoTDs inside and outside the circular trajectory. However, the period T set for UAV trajectories may invalidate the use of r cp 2 as a radius if π r cp > v max T . For such a case, the maximum valid circular trajectory radius, r max , is obtained as r max = v max T 2π . Hence, the initialized circular trajectory radii are set to r trj = min r max , r cp 2 . We set the initialized angular trajectory as θ n 2π (n−1) N −1 , ∀n. Then, with the determined values for c m trj and r trj , the initial position of UAV m in time slot n is obtained as follows: trj + r trj cos θ n , y m trj + r trj sin θ n ] T , ∀m, n. (31) Hence, we obtain the initial trajectory set Q 0 as Q 0 = {q 0 m [n], ∀m, n}. The results for the case with M = 2 and M = 3 are illustrated in Figure 2.
For the initial UAV-IoTD associations and IoTD transmit powers, an arbitrary choice is made and quickly checked for initial feasibility with regards to problem P 0 coupled with the initial circular trajectories obtained above. We set all initial powers to P max /2 and distribute the service time evenly between the IoTDs. However, this is more of an arbitrary choice that is verified for initial feasibility; many other possible combinations of initial values can find convergence.

B. SIMULATION PARAMETERS
We assume the IoTDs are uniformly and randomly distributed across a 2 × 2 km 2 area, each of them having a task to be  computed through offloading to UAV-mounted MEC servers. We assume that the IoTDs have identical energy budgets and required CPU cycles per bit, i.e., E k =Ē and C k =C. In all cases, the total bandwidth available for uplink transmission is shared equally and orthogonally between the UAVs. The UAVs are further assumed identical in terms of battery type and capacity. All the relevant parameter settings are listed in Table 1 [10], [21], [22], [23], and [27].

C. CONVERGENCE OF THE PROPOSED ALGORITHM
The first step in validating the effectiveness of the proposed scheme is to ensure the convergence of the P-BCD-based algorithm within a practical number of iterations. The penalty parameter λ 0 is initially set to a very small value, λ 0 = 10 −4 , to eliminate the unwanted influence of the initial set of binary UAV-IoTD association values on the final solution. In the outer loop of Algorithm 2, the increase in the penalty parameter with each iteration is set as c = 10, to reach large penalty effects relatively quickly. Figure 3 shows the trend of the obtained sum bits offloaded from all IoTDs with respect to the number of iterations. The curve is a result of averaging the sum bits obtained for 50 random realizations of the IoTD locations. We note that the zeroth iteration result represents the sum bits value obtained by direct calculation based on the initializations adopted. It can be seen that the proposed algorithm is quite efficient as a convergence is attained to within the prescribed accuracy = 10 −4 after only 8 iterations.

D. COMPARISON WITH BENCHMARKS AND OMA
First, to illustrate the performance gain of our proposed scheme, in Figure 4, we compare its performance to a scheme where optimization of UAV-IoTD associations with fixed circular trajectories according to the initialization scheme and full IoTD transmit powers are used. This is done versus different offloading periods T . We first observe that, as expected, the sum bits offloaded in both cases gets larger with increasing T as the IoTDs are offloaded more bits in the added time. Second, by comparing the two schemes, we observe the performance gain obtained by utilizing the proposed trajectory design. Further, this gap increases with increasing T . This is because with larger T , the optimization of UAV trajectories takes a bigger priority in order to achieve better links with the offloading NOMA IoTD clusters. In summary, joint optimization of the UAV-IoTD associations, the transmit powers, and the UAV trajectories ensures that we obtain a significantly higher value for the offloaded sum bits compared to other cases optimizing a subset of the decision variable set.
Second, to further prove the advantage of our proposed model, we do a direct comparison with its OMA counterpart. In particular, we use the TDMA version of our system. In the TDMA system, in addition to each IoTD being associated  with at most one UAV in any time slot, we further restrict each UAV to be associated with at most one IoTD in each time slot. The resulting rate expression can be formulated as We start by showing the horizontal projections of the optimized trajectories of the UAVs for the proposed NOMA scheme and the reference OMA (TDMA) scheme in Figure 5. Table 2 compares the sum bits offloaded and average UAV energy consumption of the 2 schemes quantitatively. In addition to the offloaded data size, UAV energy consumption also plays an important role as a performance metric considering the reliance on batteries. From the table, we note that our proposed NOMA scheme increases the sum bits quantity by 7.9% and reduces the average total energy consumption per UAV by 43.7%. In this scenario, and generally, the ability of all IoTDs to fully utilize the time and frequency resources to offload their data in NOMA with reduced interference due to SIC works to boost the sum bits offloaded in the system before the deadline T compared to the traditional TDMA case where the time resource cannot be shared and is optimally partitioned between the IoTDs for offloading. Further, in NOMA, the dynamic grouping/clustering of the IoTDs for offloading means that the UAVs do not need to move as much as in the case where the IoTDs are individually served at a time. This is since, for the latter case, in each time slot, each UAV would focus on the individual IoTD it is serving and would hence attempt to move closer to it to achieve better links within trajectory time and speed constraints, while for the former case, each UAV has to look at collective offloading of a number of IoTDs and does not need to move around as much to extract the most amount of offloaded bits from the IoTDs. Note that in Figure 5 the trajectory of UAV 1 for the NOMA system approximately traces a straight line between the 2 IoTDs enclosed by the initial trajectory. Further, this trajectory is almost superposed by a section of the trajectory of UAV 1 for the OMA case.
To further solidify the superiority of the NOMA case over the OMA (TDMA) counterpart, we compare the averaged obtained values for the offloaded sum bits and the energy expenditures per UAV for both the NOMA and OMA cases for 50 random realizations of the IoTD locations. The results are shown in Table 3 where a 5.8% increased offloaded sum bits for the NOMA case over the OMA case is obtained, combined with a 26.9% reduced total energy expenditure per UAV for NOMA over OMA.
Next, we compare the offloaded sum bits achieved by 3 trajectory design schemes as follows: 1) Proposed trajectory, obtained by implementing Algorithm 2; 2) Circular trajectory, obtained by implementing the proposed initialization scheme with M = 2; 3) Static UAVs, obtained by fixing UAV m at c m trj in the initialization scheme. For all three schemes, both the UAV-IoTD associations and IoTD transmit powers are optimized by Algorithm 2 with the corresponding trajectory design scheme. We further compare the results from these 3 schemes with the reference OMA scheme with joint optimization of all the decision variables (including the trajectory). The plots are shown in Figure 6. The superiority of the proposed trajectory design scheme can be observed over the others. The importance of proper trajectory design is made clear as while the NOMA case with the proposed trajectory design performs better than its OMA counterpart, the cases with preset trajectories perform worse than the OMA case with optimized trajectory. Interestingly, the case with static UAVs performs better than the one with circular trajectories; this can be attributed to the disposition of NOMA to clustering the IoTDs based on their relative positions and hence the favorability of positioning the UAVs close to their cluster centers, where the centers of the packing circles in the proposed initialization scheme serve as a fair basic way of centering the UAVs. As explained before, the gap between the different schemes increases with increasing T due to the higher impact of proper trajectory design on the obtained sum bits.

E. NUMBER OF UAVs
In this subsection, the impact of the number of UAVs on the system performance is discussed. Here, we compare the offloaded sum bits and the UAV energy expenditures for the cases M = 1 to M = 3 employed UAVs, keeping the other parameters unchanged. The results are averaged over 50 random realizations of the IoTD positions. As shown in Table 4, comparing the cases M = 1 and M = 2 UAVs, we obtain a 3.8% higher offloading volume for the case with M = 2 and a 44.6% decrease in the total energy expenditure per UAV. This means that we can obtain a significantly higher energy efficiency for offloading in the M = 2 case. This can be attributed to the need for a single deployed UAV to go through a longer way to approach as many IoTDs as possible, while in the case of 2 deployed UAVs, each UAV needs to cover a portion of the IoTDs, making it more possible to approach its served IoTDs with significantly reduced distance to cover. However, when comparing the case M = 2 with the case M = 3, we see that the increase in energy efficiency and offloaded sum bits does not necessarily scale with the number of UAVs for M > 2. This is because, for the M = 3 case, the additional UAV creates unnecessary flight energy cost. This is combined with the increasing negative impact of splitting the spectrum between the 3 UAVs, compared to the case with 2 UAVs, thus undermining the spectral efficiency of NOMA. From Table 4, we see that the M = 2 and M = 3 cases result in almost the same performance in terms of both offloaded data size and energy efficiency. Hence, adding more UAVs to the M = 2 case in our system model proves to be pointless and costly. In general, the optimal number of UAVs to adopt would depend on factors such as the network scale and the topology of the IoTDs, hence such aspects of the system model need to be carefully studied when making the choice for the number of UAVs to employ in a multiple-UAV-enabled MEC IoT system.

VI. CONCLUSION
In this paper, we proposed a multiple-UAV-enabled MEC IoT system integrating NOMA for task offloading. Given the energy constraints of the devices and the QoS requirements of the IoTDs, the joint optimization of the UAV-IoTD associations, the IoTD transmit powers, as well as the UAV trajectories was performed with the aim of maximizing the sum bits from all IoTDs within the given period for offloading. To solve the formulated highly-coupled MINLP problem, a P-BCD algorithm is developed to deal with the binary constraints and to decompose the problem into 4 subproblems. The algorithm alternately solves the four subproblems in the inner loop by applying the SCA technique and iterates in the outer loop until a suboptimal solution is obtained. Numerical results were provided to demonstrate the effectiveness of our proposed algorithm in increasing the offloaded data size from all IoTDs and improving the energy efficiency compared to several benchmark schemes and the reference OMA scheme often used in the literature. Specifically, under carefully chosen practical system operational conditions, the proposed NOMA scheme is shown to increase the size of offloaded data obtained from the OMA counterpart by 5.8% and reduce the total energy expenditure per UAV by 26.9% on average for a random placement of IoTDs across the given region, marking sizeable gains in offloading efficiency and system longevity.
There are a number of additional concerns that could be considered for future work. Namely, the effect of changing the number of IoTDs that need to be served on the system performance for the NOMA scheme compared to the OMA one. Further, for large number of IoTDs in the system, limiting the number of IoTDs associated with each UAV becomes a necessity to avoid dealing with a prohibitively complex SIC process with significant error propagation. Additionally, to keep the algorithm convergence speed from being noticeably affected by increasing the number of IoTDs, a hybrid solution scheme replacing parts of the SCA P-BCD algorithm with heuristic approaches that achieve close approximations with notably faster convergence speeds could be proposed.