I Introduction
In the future Internetofthings (IoT), the IoT devices are foreseen to be widely deployed to collect the data required by a myriad of applications. However, in certain areas where the terrestrial network coverage is unavailable due to signal blockage, spectrum scarcity or low IoT transmit power, collecting the IoT data becomes a challenging issue [1]. Although such an issue can be alleviated through densely deploying massive BSs or smallcells, the prohibitive capital expenditure (CapEx) and operating expenditure (OpEx) are unacceptable for IoT operators. In addition, the IoT data traffic is spatiotemporally unbalanced and dynamicdistributed in dedicated areas [2], which leads to frequent network congestion in the case of fixdeployed networks. To address the connection and flexibility shortages, the stateoftheart drone communication technology is leveraged to support IoT data collections [3].
Various field experiments have verified the potential of drone, i.e., unmanned aerial vehicle (UAV), to support communication services for terrestrial users. Equipped with certain communication modules, the flying drone can perform as the drone cell (DC) which collects the data from IoT devices through the user (IoT device)todrone (U2D) wireless links. Compared with the terrestrial BS, DC has distinct advantages:

Lineofsight (LoS) connection: UsertoBS (U2B) wireless links in terrestrial IoT are frequently blocked by the terrain or buildings, and such nonlineofsight (NLoS) wireless links seriously degrade the U2B communications [4]
. In contrast, the flying DCs are capable of adjusting 3D positions, which can ensure higher probability to connect terrestrial IoT devices via the highly reliable LoS links
[3]. 
Dynamic deployment: Compared with terrestrial BSs built on fixed locations, DCs can change their deployments according to the spatiotemporal traffic variations, and assigned to different controllers or IoT devices on demands [5].

Fullycontrolled mobility: The flying traces and behaviors of any DC are fully controlled by the network providers, which empowers DCs with the dynamic deployment feature and facilitates interDC collaborations [6].
Leveraging the advantages of DC, how to control the 3D placements of DCs to improve the network performance is of great importance yet very challenging due to the complicated 3D mobility of DC and dynamic DCtoground (D2G) link quality. In the literature, there have been several studies on the 3D placements of single/multiple DCs to support terrestrial IoT, which can be classified as two categories, i.e., static DC deployment and DC trajectory design. The static DC deployment investigates the optimal hovering positions of DCs to maximize the communication performance of associated users. BorYaliniz
et al. proposed a bisection search based algorithm to deploy single DC [7]. Kalantari et al. optimized the multiple DCs static deployment through the swarm intelligence based algorithm [8] [9]. To maximize the IoT information collection gain, Mozaffari et al. leveraged the clustering based method to design optimal hovering positions for multiple DCs [10]. However, the static deployment can hardly guarantee the fairness among users. Particularly, the IoT devices located at the edge of the DC’s radio coverage can suffer severer pathloss compared with devices located at the center of the radio coverage, which leads to unbalanced communication performance.To promote the fairness among IoT devices, DC trajectory design schemes are proposed in some recent works, in which DCs periodically fly over associated users via optimized trajectories, and communicate with dedicated users in scheduled slots [12][16]. Following predetermined DC trajectories, Li et al. investigated the packets transmission relayed by multiple DCs from the sensor to BS [11]. Zeng et al. proposed a general framework for trajectory and communication optimization of a energyefficient pointtopoint UAVground communication system [12]. In [13], Wu et al.. formulated the nonconvex DC trajectory design problem in which the user association, DC horizontal trajectory design and DC transmitting power control are jointly optimized; the delay constraint is further considered in the extension work focusing on delaysensitive services [14]. However, there are some strong assumptions and model simplifications in current DC trajectory design studies. For instance, in current works, the flying altitudes of all DCs are usually set to be a constant value to simplify the optimization process (i.e., only 2D trajectory is optimized), and the adopted D2G channel model is the classic Friis model which fails to characterize the distinct features of DCs [13] [14]. Besides, most existing trajectory design works idealize or ignore the drone DCtoBS (D2B) link quality constraints, which however is essential for droneassisted networks.
To address the aforementioned issues, in this paper, we study the multiple DCs 3D trajectory design to facilitate the IoT data collection. We leverage the stateoftheart U2D [15] and D2B [16] pathloss models, and the DC altitudes are taken as optimization variables to further improve the network performance. We formulate the 3D multiDC trajectory design as a mixed integer nonlinear programming (MINLP) problem, in which the summation of average U2D pathloss is minimized under the constraints of D2B link qualities, user fairness and horizontal/vertical speed limitations. Given the constant transmit power, bandwidth and interferencelevel of IoT devices, the lower average U2D pathloss indicates higher network throughput. As the formulated MINLP problem is very difficult to be solved directly, we decouple the original MINLP problem into multiple quasiconvex or linear subproblems. Each subproblem optimizes the user association, user scheduling, DC horizontal trajectories and flying altitudes, respectively. By adopting the block coordinate descent (BCD) method, we then propose a 3D multiDC trajectory design algorithm in which those subproblems are iteratively optimized to achieve the optimal result of the MINLP problem. According to the simulation results, our algorithm can reduce the average D2U pathloss by , and lower the standard deviation by 56% when compared with the static DC deployment scheme.
The remainder of this paper is organized as follows. In Section II, we introduce the system model. Then, we formulate the trajectory design problem for multiple DCs in Section III. In Section IV, we propose the 3D multiDC trajectory design algorithm to solve the problem. Numerical results are shown in Section V, and the conclusion are given in Section VI.
Ii System Model
The overview of DC enabled IoT data collection is shown in Fig. 1, where multiple DCs are released to relay data uploading from IoT devices to one BS. We define the DC users as the IoT devices whose reliable U2B links are compromised due to blockage or resource shortage. The BS periodically detects the channelstateinformation (CSI) of each U2B pair and assign compromised IoT devices into the DC user set . Given and the available DCs set , the BS implements the trajectory design algorithm to determine the optimal trajectory for each DC. DC users are considered as identical devices which transmit data to the BS with identical data collection period, bandwidth, and transmit power. Since the proposed trajectory design algorithm can schedule the transmit time of all users associated to one DC to prevent transmitinterval overlapping, and we assume each DC can distinguish its associated users’ signal efficiently from others, the interuser interference issue is prevented in our system. The cardinalities and represent the number of DC users and available DCs, respectively.
Iia U2D and D2B Channel Models
Both U2D and D2B wireless links are modelled based on the recent D2G channel research [15] [16]. According to [15], we can express the U2D pathloss as
(1)  
where is the U2D horizontal distance, and is the DC flying altitude. and are carrier frequency and speed of light in Hz and m/s, respectively. and are environmentdependent LoS and NLoS pathloss offsets, respectively. is the U2D LoS probability, which is expressed as
(2) 
where and are environmentdependent parameters. To prevent the interference to U2B communications, as well as to provide additional spectrum resource for DC users, is expected to be different from the traditional IoT band, such as WiFi or NBIoT [17] bands. The average D2B pathloss is calculated by [16]
(3) 
where and denote the D2B horizontal distance and the DC to BS antenna vertical angle, respectively. , , , , and are all environmentdependent parameters. Since the 850 MHz band is used for LTE transmissions [16], (3) contains no parameter representing carrier frequency.
IiB DC Trajectory Model
For one DC , we consider that it serves the associated user set through a time division multiple access (TDMA) mode. Within one trajectory period , DC flies over all its associated users and sequentially severs them according to the scheduling result. is set to the same value as the data collection period of users, in order that each user can transmit the collected data once within one period. By discretizing into equaltime slots, trajectory of DC within one period can be modeled as an
length sequence composed by threedimensional vectors:
(4) 
where , and denote the 3D coordinates of DC at slot . is the set of slots in one period. In this work, we define when user is associated to DC for , otherwise
. The U2D communication scheduling is denoted by the binary variable
for . If user is severed by DC in slot , is set as ; otherwise, .Several trajectory constraints are considered in our work: 1) The trajectory of each DC has to be a 3D closed curve since the DC must return to the start point by the end of each period for a new period. 2) Since the power consumption models for horizontal and vertical moving of DC are different, it is reasonable to set the maximum horizontal speed and maximal vertical speed , respectively. 3) In any slot , one DC can serve at most one user ; in all slots, one user can only be associated to one DC. 4) The slots amount scheduled to one within one period cannot be smaller than a predefined threshold , which indicates the minimal time to complete each U2B data transmission.
Iii Problem Formulation
We define as the coordinate of , as the 2D projection of DC on XY plane, and as the horizontal distance between and . Define the elevation angle between DC and BS antenna as in degree. is the Euclidean norm.
By substituting and into (1) and (3) respectively, we can calculate the U2D pathloss between DC and user in slot by
(5)  
as well as the D2B pathloss between the BS and DC by
(6)  
Assuming all users maintain their transmission power in each period, the effectiveness of the data collection for each U2D link is determined by the U2D pathloss. Therefore, we aim to minimize the total pathloss suffered by each DC in the trajectory design. Define the matrices , and . The 3D multiDC trajectory design problem can be formulated as
(7a)  
(7b)  
(7c)  
(7d)  
(7e)  
(7f)  
(7g)  
(7h)  
(7i)  
(7j)  
(7k) 
In (7), is the maximum number of users for one , and is the D2B pathloss threshold. (7c)(7g) correspond to user association constraint 3) and 4) in section II. (7h)(7j) represent the DC trajectory constraints 1) and 2) in section II. (7k) is the D2B pathloss constraint.
Due to the quadratic and exponential terms in (7a) and (7k), and the binary variable and , the 3D multiDC trajectory design problem turns to be an MINLP problem [7]. Besides, the objective function (7a) and constraint (7n) are both nonconvex for DC trajectory , which increases the difficulty to solve the problem.
Iv 3D MultiDC Trajectory Design Algorithm
The nonconvex problem (7) can be transformed into solvable forms (e.g. quasiconvex or LP) by assuming parts of the decision variables as constants. Then, we can decouple the MINLP problem into multiple subproblems that are solvable for parts of the decision variables, and iteratively solve them by applying the BCD method. Specifically, we divide the decision variable in problem (7) into four blocks (i.e. , , and ), and formulate following subproblems.
Given the predefined trajectories of multiple DCs (constant , and ), the user association subproblem can be written as an integer linear programming (ILP) problem:
(8)  
Problem (8) can be efficiently solved through the branch and bound method, which is well supported by Gurobi solvers.
Based on the optimized and predefined trajectory, the U2D communication scheduling subproblem is also transformed to an ILP problem:
(9)  
Same as problem (8), the Gurobi optimizer can efficiently solve problem (9).
According to (6), for given , and , (7a) is still nonconvex for . However, by further decoupling into element optimization variable , (7a) is both quasiconvex and nondecreasing function to U2D horizontal distance , when other keep constants. Since (7a) can reach its minimal value when for any within its available range, the objective of minimizing (7a) equals minimizing , which is a quadratic convex function for . On the other hand, [18] proves that the available range of for constraint (7k) forms a convex set in XY plane. Therefore, the subproblem to find optimal can be formulated as
(10)  
whose optimal can be calculated by solving the following quadratic convex optimization problem:
(11)  
For any given , and , the subproblem optimizing is also nonconvex to due to the form and a form terms. Same as preceding decoupling process of , can also be decoupled into element optimization variable , and the subproblem to optimize is
(12)  
Problem (12) can only be proved quasiconvex to when other keep constant. However, since is the summation of one nonincreasing function of and one nondecreasing function of , we can argue that problem (12)’s objective function has only one global minimum. To achieve the optimal which minimizes problem (12), we can leverage the NewtonRaphson method by iteratively calculating the following function:
(13) 
where the iteration stops when and .
The subproblems of problem (7) can all be optimized respectively with other optimization variables keeping constant. Therefore, problem (7) can be solved through iteratively optimizing those subproblems until the results converge, which yields to the classic BCD method. Based on the BCD method, we propose the 3D multiDC trajectory algorithm as shown in Algorithm 1. , , denote the multiple DCs’ user association, U2D communication scheduling and DC trajectories calculated after iteration , respectively. Since the global optimal results of each subproblem can be achieved accurately, the proposed algorithm ensures convergence [13] [19].
V Simulation Results
We implement the proposed algorithm in simulations by using Gurobi solver [20]. The parameters of both U2D and D2B pathloss models are set as suburban scenario [16]. The carrier frequency for U2D communication is set as , which is widely supported by commercial drone products and prevents interference to the cellular system. D2B communications use the LTE band [16]. As the initial trajectories, DCs are uniformly deployed over the radio coverage area of terrestrial BS with the same flying altitude . Without loss of generality, we treat one slot as the minimal time unit to calculate related variables including , and . According to the general specifications of commercial drones whose maximal horizontal and ascent/descent speeds are and respectively, the approximate value of is around . Table. I summarizes the detail simulation parameters.
Parameter Name  Value 

BS radio coverage radius  
User(IoT device) number  
U2D parameters  (0.1,21,4.88,0.43) 
D2B parameters  (3.04,23.29,3.61,4.14,20.7) 
Carrier frequencies (U2D, D2B)  
Duration of one period  
D2B pathloss constraint  
Minimal number of slots  
Maximal horizontal speed  
Maximal vertical speed  
Trajectory difference  for each slot 
Fig. 2 shows an example of trajectory design result in which five DCs relay data from 20 IoT users. The closed curves with different markers in Figs. 2(a) and 2(b) denote the trajectories of different DCs, while the squares on the XY plane represent the users. Users are associated to their corresponding DCs with same colors. Fig. 2(c) shows the variations of flying altitude within one period, where the lower bound of all trajectories near is constrained by the D2B link quality constraint. As shown in Figs. 2(a) and 2(b), for each DC, the optimized trajectory can fly over all its associated user and form a closed trajectory in 3D space. In Fig. 2, the time interval between two adjacent dots on the same trajectory is exactly one slot.
In Fig. 2(b), we can see that more than dots on each trajectory are overlapped above the associated users positions. By joint analyzing Figs. 2(b) and 2(c), we can justify that DCs are prone to hovering above the associated users, while spending less slots for the travel process between two hovering positions. Such a “hovering effect” in the optimized trajectories can indicates the effectiveness of our proposed algorithm in minimizing average U2D pathloss. Because with given flying altitude, the average U2D pathloss can be minimized when U2D horizontal distance is zero, i.e., hovering above the user.
Since most existing trajectory planning results are based on different D2G channel models with fixedheight (i.e. [13]
), we choose one static DC deployment algorithm, the perdrone iterated particle swarm optimization (DIPSO) algorithm
[18], as the benchmark to highlight the efficiency of our proposed algorithm. The average U2D pathloss performance achieved by the proposed 3D multiDC trajectory design algorithm, as well as the static DC deployment scheme are compared in Fig. 3 and Fig. 4. In both figures, the blue curves represent static DC deployment while the red curves represent DC trajectory design.From Fig. 3, we can see that the U2D pathloss performance of both algorithms are improved as the available DC number increases. However, the average U2D pathloss of our trajectory design solution remains less than that of the DIPSO solution by . The CDF curves in Fig. 4 further confirm that the trajectory design algorithm can maintain the U2D pathloss below given thresholds with higher probability when compared with static DC deployment.
The user fairness promotion provided by the DC trajectory design algorithm is indicated by the error bars in Fig. 4 and the U2D pathloss standard deviation comparison in Table II. and are U2D pathloss standard deviations for DC trajectory design and static DC deployment, respectively. In Fig. 4, we can observe that the average U2D pathloss of static deployment ranges from to , which is four times as much as that of the trajectory design algorithm. From Table II, we can conclude that the trajectory design can reduce the U2D pathloss standard deviation by on average compared with static deployment.
DC number  

Vi Conclusion
In this paper, we have investigated the 3D multiDC trajectory design for efficient IoT data collection. An MINLP problem has been formulated to minimize the summation of average U2D pathloss. To solve the MINLP problem, we have decoupled it and formed multiple subproblems in which the user association, U2D communication scheduling, horizontal trajectories, and flying altitudes are optimized, respectively. Leveraging the BCD method, we have devised the 3D multiDC trajectory design algorithm to solve the MINLP problem by solving subproblems iteratively. Simulation results have shown that the proposed DC trajectory design algorithm can achieve average U2D pathloss reduction, and promote pathloss standard deviation by more than when compared with the static DC deployment. In future works, we will analyze the impacts of initial deployments, horizontal and vertical flying speeds, as well as interDC safe distance on the trajectory design, and investigate the communication and computation resources allocation of multiple DCs with the optimized trajectories.
Acknowledgment
This work is supported by the National Natural Science Foundation of China under Project 91638204 and the Natural Sciences and Engineering Research Council (NSERC), Canada.
References
 [1] N. Cheng, W. Xu, W. Shi, Y. Zhou, N. Lu, H. Zhou, and X. Shen, “Airground integrated mobile edge networks: Architecture, challenges, and opportunities,” IEEE Commun. Mag., vol. 56, no. 8, pp. 26–32, 2018.
 [2] M. Tao, K. Ota, and M. Dong, “Locating compromised data sources in IoTenabled smart cities: A greatalternativeregionbased approach,” IEEE Trans Ind. Informat., vol. 14, no. 6, pp. 2579–2587, 2018.
 [3] M. Mozaffari, W. Saad, M. Bennis, Y.H. Nam, and M. Debbah, “A tutorial on UAVs for wireless networks: applications, challenges, and open problems,” 2018. [Online]. Available: https://arxiv.org/abs/1803.00680.
 [4] J. Liu, Y. Shi, Z. M. Fadlullah, and N. Kato, “Spaceairground integrated network: A survey,” IEEE Commun. Surveys Tuts., early access, 2018.
 [5] Z. Xue, J. Wang, G. Ding, H. Zhou, and Q. Wu, “Maximization of data dissemination in UAVsupported internet of things,” IEEE Wireless Commun. Lett., early access, 2018.
 [6] S. Zhang, W. Quan, J. Li, W. Shi, P. Yang, and X. Shen, “Airground integrated vehicular network slicing with content pushing and caching,” IEEE J. Sel. Areas Commun., early access, 2018.
 [7] R. I. BorYaliniz, A. ElKeyi, and H. Yanikomeroglu, “Efficient 3D placement of an aerial base station in next generation cellular networks,” in Proc. IEEE ICC, Kuala Lumpur, Malaysia, May 2016, pp. 1–5.
 [8] E. Kalantari, H. Yanikomeroglu, and A. Yongacoglu, “On the number and 3D placement of drone base stations in wireless cellular networks,” in Proc. IEEE VTCFall, Montreal, QC, Sept 2016, pp. 1–6.
 [9] E. Kalantari, M. Z. Shakir, H. Yanikomeroglu, and A. Yongacoglu, “Backhaulaware robust 3D drone placement in 5g+ wireless networks,” in Proc. IEEE ICC Workshops, Paris, France, May 2017, pp. 109–114.
 [10] M. Mozaffari, W. Saad, M. Bennis, and M. Debbah, “Mobile unmanned aerial vehicles (UAVs) for energyefficient internet of things communications,” IEEE Trans. Wireless Commun., vol. 16, no. 11, pp. 7574–7589, 2017.
 [11] K. Li, W. Ni, X. Wang, R. P. Liu, S. S. Kanhere, and S. Jha, “Energyefficient cooperative relaying for unmanned aerial vehicles,” IEEE Trans. Mobile Comput., vol. 15, no. 6, pp. 1377–1386, 2016.
 [12] Y. Zeng and R. Zhang, “Energyefficient UAV communication with trajectory optimization,” IEEE Trans. Wireless Commun., vol. 16, no. 6, pp. 3747–3760, 2017.
 [13] Q. Wu, Y. Zeng, and R. Zhang, “Joint trajectory and communication design for multiUAV enabled wireless networks,” IEEE Trans. Wireless Commun., vol. 17, no. 3, pp. 2109–2121, 2018.
 [14] Q. Wu and R. Zhang, “Common throughput maximization in UAVenabled ofdma systems with delay consideration,” 2018. [Online]. Available: https://arxiv.org/abs/1801.00444
 [15] A. AlHourani, S. Kandeepan, and S. Lardner, “Optimal LAP altitude for maximum coverage,” IEEE Wireless Commun. Lett., vol. 3, no. 6, pp. 569–572, 2014.
 [16] A. AlHourani and K. Gomez, “Modeling cellulartoUAV pathloss for suburban environments,” IEEE Wireless Commun. Lett., vol. 7, no. 1, pp. 82–85, 2018.
 [17] M. E. Soussi, P. Zand, F. Pasveer, and G. Dolmans, “Evaluating the performance of eMTC and NBIoT for smart city applications,” in Proc. IEEE ICC, Kansas City, MO, May 2018, pp. 1–7.
 [18] W. Shi, J. Li, W. Xu, H. Zhou, N. Zhang, S. Zhang, and X. Shen, “Multiple dronecell deployment analyses and optimization in drone assisted radio access networks,” IEEE Access, vol. 6, pp. 12 518–12 529, 2018.
 [19] D. P. Bertsekas, Nonlinear programming. Belmont, MA, USA: Athena scientific, 1999.
 [20] “Gurobi optimizer 8.1.” [Online]. Available: http://www.gurobi.com