Deployment Strategies of Multiple Aerial BSs for User Coverage and Power Efficiency Maximization

1 Deployme nt Strate gies of Multiple Aerial BSs for User Co v erage and Po wer Ef ficienc y Maximizat ion Jingcong Sun, Studen t Member , IEEE, Christos Masouros , Senior Member , IEEE Abstract —Unmanned aerial ve hicle (U A V) based a erial base stations (BSs) can pro vide rapid communication services to ground users and are thu s promising f or future communication systems. In this paper , we consider a scenario where no functional terrestrial BSs are a vailable and the aim is deploying multip l e aerial BSs to cov er a maximum numb er o f users within a certain target area. T o this end, we first p ropose a naiv e successive deployment method, which con verts the non-con vex constraints in the i n vo lved optimization into a combination of l inear constraints through geometrical relaxation. T h en we inv estigate a deployment method based on K -means clustering. Th e method divides the target area into K conv ex subar eas, where within each subarea, a mixed integer non-linear problem (MINLP) is solved. An iterativ e power efficient technique is further proposed to impro ve cov erage probability with reduced po wer . Finally , we pro pose a r obust technique fo r compensating the loss of cov erage probability in the e xistence o f inaccurate user location i nform ation (ULI). Our simulation results show that, the proposed techni ques achiev e an u p to 30% higher cov erage probability when users are n ot distributed unif ormly . In addition, th e pro posed simultaneous deployment techn iques, especially the one using i t erative algo- rithm impro ve power -efficiency by up to 15% compared to the benchmark circle packing theory . Index T erms —Unmanned aerial vehicles, user coverage, air - to-ground comm unication, clu stering algorithm I . I N T RO D U C T I O N L O W -AL TITUDE unmann ed aerial vehicles ( U A Vs) have been increasingly ap pealing to future wireless commun ication systems. U A Vs which are cost-effectiv e, interoper able, flexible and likely to have a higher probab ility of line - of-sight (L oS) ch annels ar e p romising in a various scenarios for bo th ci vilian and military use [1]–[3]. U A V -based commun ication ap plications in volve three m ain categories, namely relaying, information dissemination/d a ta co llection and ubiquitous c overage [4]. U A Vs serving as relaying nodes ar e studied to provide reliable w ir eless com munication between distant users with blocked direct lin ks while increasing the system throug hput [5]–[ 8]. Optimizin g the flying trajector y is o f particular interest when UA Vs are dispatched to d isseminate inform ation [9], [10]. L ast but not least, as U A Vs can be q uickly dep loyed, U A V -based aerial base stations (BSs) are attracting in c r easing interest as a mean s to p rovide fast wireless services to gro und users. For instance, UA Vs can b e dep loyed to ease th e burden of terrestrial ba se station s in extremely crowded areas by offloading u sers fro m gro und cells when specific rate or distance requ ir ement is satisfied [ 11], [12]. Moreover , fast J.Sun and C.Masouros are wit h the Department of E lectronic a nd Electrical Engineeri ng, Uni versity Coll ege Londo n, London WC1E 7JE, U. K. (e-mail: uceej su@ucl.ac.uk; c hris.masouros@ie ee.or g) service recovery or supply can be offered by such flying BSs when fixed communication infrastructure s are damaged. W ith the rising in terest in the above areas, the c hallenges in the p ractical use of aerial BSs are becoming pertinen t. When the aerial BSs are utilized for em ergency commu nications such as search -and-r e scu e, the pr iority is finding the optim al locations of UA Vs so that a maximum number of u sers can be covered. Meanwhile, since built-in batteries are used for supplying power in most cases, limited on - board energy is another factor th a t co n strains the end urance of aerial BSs [4] , [13]–[15]. It has been proven that p rolong ed op eration time can be achieved by redu c in g the transmit p ower o f aerial BSs when qua lity-of-serv ice (QoS) requir ements a r e met [1 6], [17]. The ae r ial BS coverage prob lem is first stud ied in [ 1 8], which gave an a ir-to-gr ound (A TG) chan nel mo del used to find the optimal altitude of U A Vs that can lead to maximum coverage area on the gro und. Instead of ju st m aximizing the coverage ar ea, recent re search h as beco me incr easingly focused on algorithms trying to cover the maximum nu m ber of user s [1 7], [19]–[21]. Specifically , [17] solved a 3- D circ le placement p roblem by formu lating it as a mixed integer non-lin e ar problem (MINLP), while [19] made a furth er step by conside r ing explicit QoS co nstraints. I n [21], optim al location of an aer ial BS is obtain ed thr ough an exhau sti ve search in pr edefined gir d s. However , all th e work s men tioned above consider ed on ly the c ase of a sing le flyin g BS which limits their use. In most real situations, it is necessary to deploy mu ltiple U A Vs at the same time to cover a majority of users in a specific target region. The work in [22] exten d s the number of utilized U A Vs to tw o with a carefu l co nsideration of inter-cell inter f erence (ICI). Mozaffari et al. [2 3] prop o sed a circle packing technique so that the total ar ea covered by multiple aerial BSs is maximized , howev er the method did not consider user d istributions. A 100% user coverage p robab ility is shown thr ough a spiral algo rithm in [24]. H owever , th e study ignores the effect of I CI, and th e interferen ce issue needs to be tackled with overlaid techn iques. In this pap er , we study th e efficient de ployment o f multiple U A Vs so the max imum user coverage probability is ach iev ed, where we define th e coverage probab ility as the r atio of number of covered users to the total numb er of u sers within a specific target area. F ollowing [23], [24], we assume that the locations o f users are known with the help of hig h-accur a cy GPS systems and each aerial BS has enoug h cap a city to supply all the u sers it covers. W e con sider a scenario where multiple UA Vs are deployed in a target area with out g r ound BSs’ coverage. This is relevant in ru r al ar ea coverage in 2 cases wher e terrestrial BS s are absent, and in natural disaster scenarios where ter restrial coverage is disabled . Rotary-wing U A Vs which h av e the ab ility to move in arbitrary direction as well as h old still in the air are assumed as the carrier for aerial BSs [ 4]. O u r aim is to stud y the efficient placem ent of m ultiple aerial BSs in order to obtain a m aximum user coverage pro bability wh ile completely a voiding the influence of ICI. The UA V placemen t problem is mo delled as a circle placement pr oblem and no coverage overlap is allowed so that ICI between aerial BSs is intrinsically av oided. Our simu lation results demon strate that the pro posed circle placement metho ds achieve higher user coverage p robab ility than the benc hmark circle p acking the ory (CPT) [23]. Moreover , the increased coverage pro bability is achieved with significantly reduced tran sm it power in certain scenarios. The existence of inaccurate user lo cation inf o rmation (ULI) is also co nsidered, and clearly increased ro bustness against inaccurate ULI is obtained when proposed robust techn iq ue is app lied. For clarity , W e summarize the ma in contributions of this paper as follows • Geome trical Relaxation : we pr o pose a ge o metrical r elax- ation scheme where the optimal locations of U A Vs are obtained in a ’ step- by-step’ fashion, in which the next U A V is always dep loyed in a p osition such that it covers the most nu m ber of rema in ing users in the target area until there is n o spac e f or acco mmoda tin g more UA Vs. The formulation of such a problem includes a increased number of non- conv ex constraints to av oid interfer e nce between a ny two aer ial BSs. Th e n on-co nvex constraints are add ressed with a simp le geom etrical relaxation which conv erts each non -conve x con straint into four linear con- straints th a t can be easily solved. • K -means Deploymen t: a more efficient te c h nique which can be easily applied to any size o f target area is proposed with the help of K -mea n s clu stering algo rithm [25], where the best lo c ations of a e rial BSs a r e fou nd within se veral subareas which are convex r egions. • Power Efficient K -mean s Deploymen t: we then propose an iterative alg orithm to further impr ove the u ser cov- erage probability while dr astically r educing the r equired transmit power of aerial BSs. • Robust Deploymen t with imperfect ULI : the coverage perfor mance in the existence of im perfect ULI is also considered , and a rob ust tech n ique to compensate the p e r- forman ce loss with minimum tr ansmit power is proposed accordin g ly by shif ting the loca tion of aer ial BSs within a b ound e d circular region before incr easing the r adii of correspo n ding coverage areas. • Complexity Analy sis: we d eriv e th e computational com- plexity of the proposed techniq ues analytically in terms of the floating-p oint opera tio ns requ ired. I t is shown mathe- matically that the comp lexity of the K - means deploymen t scales linearly with the numb er of UA Vs a nd the num b er of users, wh ile an add itional q uadratic scaling with the number o f UA Vs is observed for the robust techniqu es. The remainde r of this paper is organiz e d as fo llows. Sec- tion II introdu c es system m o del. Th e successi ve dep loyment method based o n geometrical relaxation is p roposed in Sectio n III. Section I V and Section V propose two simultaneo us de- ployment metho ds based on K -means c lu stering, and Sectio n VI consider s the situation o f inaccurate ULI. Compu tational complexity of the proposed techniques are analyzed in Section VII and numeric a l resu lts are shown in Section VIII. Finally , the pap er is concluded in Section IX. I I . S Y S T E M M O D E L W e consider a squar e geo g raphical target ar ea with side length L s containing a set of low-mobility users denoted b y M as shown in Fig. 1. W e assume K multiple aerial BSs are deployed within the region in or d er to provide wir eless commun ication to as many gr ound users as p ossible. No te that, as U A Vs can move fr eely , such d eployment of aerial BSs can be done regularly in o rder to accommodate any c hanges in the user po sitions. W e no te that trajector y design is ou t of the scope of this paper and we will only fo cus o n each snapsho t of users within the ar ea. A. P ath lo ss mod e l There are a few different models char acterizing air-to- groun d (AtG) links. Several paper s simply u se LoS channels as the AtG ch annel model due to the fact tha t the U A V - groun d commu nication is more likely to be dominated by LoS links compa r ed to terrestrial com m unication systems, especially in areas with o ut a large number of high-r iser s such as suburban and rural areas [10], [24]. Her e , we utilize a mo re strict and gener a l ch annel m odel as in [18] which con siders both the ef fect of LoS link s and non- lin e-of-sigh t (NLoS) links. I f we de n ote the location o f user i in the set M as ( x i , y i ) , the horizo ntal location and the altitu de of the k -th U A V as ( x ck , y ck ) a n d h k , k = 1 , 2 , ...K respectively , th en the ground d istance between user i and U A V k is r ik = q ( x i − x ck ) 2 + ( y i − y ck ) 2 . Acco rdingly , the probability of having a LoS link is: P LoS = 1 1 + a exp( − b ( 180 π θ ik − a )) (1) where a and b ar e parameters depending o n the en vironmen t and θ ik = tan − 1  h k r ik  is the elev a tio n ang le in radians. Note that the pr obability of having a NL oS link is repre sen ted by P NLoS =1 − P LoS . In practice, path loss may b e affected by specific topolo gy , groun d morpho logy , shad owing effect as well as rap id flu c- tuations over short time duration s. Ho we ver , following [17]– [19], we con sider the m e a n effect of path lo ss instead of its small-scale fluctuations in this pap e r . As ind icated in [1 8], th e mean p ath loss can be modeled a s free sp ace p ath loss plus a n excessi ve p a th loss wh ic h depend s o n the propagation grou p. 3 Fig. 1. System m odel Therefo re, the path loss mo del for LoS an d NLoS links in d B are respec tively shown as PL LoS = 20 lo g  4 π f c d ik c  + η LoS PL NLoS = 20 lo g  4 π f c d ik c  + η NLoS (2) where f c is the car rier frequency of the system, c is the light speed, η LoS and η NLoS are the excessi ve loss for LoS and NLoS lin ks respectiv ely . In addition, d ik represents the Euclidean distance b e twe en user i and the k -th aer ial BS, which is given by d ik = p r ik 2 + h k 2 . T he weighted path loss between user i and the k -th aerial BS, wh ich is a fu n ction o f both h k and r ik is then expressed as PL ( h k , r ik ) = P LoS PL LoS + P NLoS PL NLoS (3) By sub stituting (1) an d (2) in to (3) , the above equation yield s PL ( h k , r ik ) = A 1 + a exp( − b ( 180 π tan( h k r ik ) − a )) + 10log  h k 2 + r ik 2  + B (4) where A = η LoS − η NLoS and B = 20 log( 4 π c ) + 20 log( f c ) + η NLoS . Alter n ativ ely , equ ation (4) can b e wr itten as PL ( θ ik , r ik ) = A 1 + a exp( − b ( 180 π θ ik − a )) + 20log( r ik cos( θ ik ) ) + B (5) In this paper , the service thr eshold is defined in terms of the received power of user i , wh ic h is deno ted by P i r . I f the transmit power of the k -th aerial BS is P k t , we ha ve P i r = P k t − PL ( h k , r ik ) (6) When the received power o f user i with regard to the k - th aerial BS is larger than or equal to th e th reshold value P min , u ser i is cov ered by the k -th aerial BS. When transm it power P k t is gi ven, this is equiv alent to saying that, user i is covered b y the k -th aerial BS when the propa gation p ath loss is smaller than or equal to a thresh o ld value γ th . Sin c e the mean effect of path loss is co nsidered in this paper, the coverage area of e ach aerial BS is a circle area with rad ius R k , wh ere R k = r ik | PL( h k ,r ik )= γ th . It has b e en shown in [17] that, the maximum coverage radius is alw a ys associated with a fixed elev ation ang le θ opt , which o n ly depends on the environmen t, and is gi ven by θ opt = tan − 1 ( h ∗ k R ∗ k ) (7) where R ∗ k is the maximum coverage radius we can obtain and h ∗ k is the correspo nding optimal altitud e o f the k - th a e rial BS that giv es R ∗ k . Based on (5) and (7), the threshold p ath loss can be written as γ th = A 1 + a exp( − b ( 180 π θ opt − a )) + 20 lo g( R ∗ k cos( θ opt ) ) + B (8) from which we can solve R ∗ k when a specific γ th is giv en. In this paper, we co nsider urban en vironmen t whose optimal elev ation angle θ opt equals 4 2 . 44 ◦ [17]. No te that when multiple U A Vs are deployed , th e interferen ce effect needs to be addressed. It is clear that ICI can be intrinsically av o ided when there is no o verlap between coverage ar e a s of aerial BSs. B. User Distribution W itho ut loss of ge nerality , in this paper, we assume a random user distribution which is o b tained throug h a spatial point p rocess (SPP). W e assum e thr ee types of SPP , namely homog eneous Poisson process, inho mogene ous Poisson pro - cess an d Poisson cluster proc e ss [26], [27]. The three SPP models are able to describ e a m ajority of user distributions in real scenar ios. Let D denote a bound ed set, X ( D ) denote a counting measure of D which calculates the random number of points in D , and µ ( D ) is a mean measur e of D , giving the expected numb er of points. 1) Homogeneous P oiss on pr ocess: In homo geneou s Pois- son pr ocess, all user points are u niform ly a n d ind epende n tly distributed within the target area W . The point den sity equals to a constant λ s , which describes the average n u mber of user points generated in a unit area. Therefore, f o r any user ( x i , y i ) , we h av e P (( x i , y i ) ∈ S ) = S W (9) for any subarea S of th e target area W . Note that th e num b er of user po ints gene r ated fo llows Poisson distribution with µ ( D ) = λ s · W , which is X ( D ) ∼ P oisson( λ s · W ) . 2) Inho mogeneous P oisson pr ocess: In h omog e neous Po is- son p rocess is a more general SPP mo del wh ich in troduc e s inhomo geneity . The constant point den sity λ s is replaced b y an intensity f unction λ ( x, y ) , which varies with locatio ns in the target are a. An exam p le intensity function can be λ ( x, y ) = c ( x 2 + y 2 ) (10) where c is a constant. Then we hav e µ ( D ) = E { X ( D ) } = Z D λ ( x, y ) dxdy (11) where E { . } is th e expectation op erator . T h e co rrespond - ing number o f generated user po ints is thus X ( D ) ∼ Poisson( µ ( D )) , with µ ( D ) obtained fr om (1 1). 4 3) P oisson cluster pr ocess: Users often gather around points of interest in real scenarios, in which case their dis- tributions inv olves clu stering. In order to describe this kind of u ser distribution, a Poisson cluster SPP is utilized [27]. Firstly , a set of ’p arent’ points S p is genera ted fo llowing homog eneous Poisson pr ocess with constan t po int density λ p . Then fo r each c ∈ S p , ’ch ild ren’ points are indep endently generated fo llowing Poisson process with inten sity function λ c ( x, y ) . Note that ’children ’ poin ts are distributed in c ircles around c o rrespon ding ’par ent’ p oints to form clusters. I I I . P RO P O S E D S U C C E S S I V E D E P L OY M E N T M E T H O D W I T H G E O M E T R I C A L R E L A X A T I O N ( S D - G R ) In this section , we p ropose a method based on succe s- si ve circle placement to find the optima l location s of a erial BSs such that a max imum user coverage probab ility can be achieved. Following [23], [2 4], a n d as shown in Fig. 2(a), we assume that a ll UA Vs h av e the same coverage radius R and are thus deployed in the same a ltitude H when a specific p ath loss req uirement is gi ven, that is h k = H , k = 1 , 2 , ..., K (12) R k = R, k = 1 , 2 , ..., K H = R tan( θ opt ) with R calculated from (8). T herefor e, placing multiple UA Vs is equ iv alent to p lac ing multiple circles in the horizo ntal plane such that the numb er of enclosed user points is max imized. U A Vs are placed in a suc c e ssi ve meth od, where at eac h step the p lacement of the aerial BS aims to cover the max imum number of rem a in ing users in the target area while e n suring that there is no overlap in coverage are as with all previously deployed BSs. The lo cation of the first UA V can be fo und by utilizing the method proposed in [17]. Let C 1 denote the coverage area of the first U A V a n d Boolean variable u i ∈ { 0 , 1 } , i ∈ M denotes the status o f user i such that the user is enclosed by C 1 when u i = 1 and is not covered by the first UA V when u i = 0 . Th en the circle p lac ement problem is formu late d as maximize x c 1 ,y c 1 ,u i X i ∈M u i (13) sub ject to ( x i − x c 1 ) 2 + ( y i − y c 1 ) 2 ≤ R 2 + M (1 − u i ) , ∀ i ∈ M u i ∈ { 0 , 1 } , ∀ i ∈ M where ( x c 1 , y c 1 ) is the hor iz o ntal location of the first UA V , i.e . the center of the coverage re gion, and M is a large constant that can be any value larger than the square of the largest distance between any two poin ts in the target area. It can b e seen th at the first constraint of (13) red uces to ( x i − x c 1 ) 2 + ( y i − y c 1 ) 2 ≤ R 2 , ∀ i ∈ M (14) when u i = 1 which is equivalent to saying that user i is covered by the first UA V , an d the o bjective function of (13 ) is increased by 1 correspon dingly . In addition, the inequality of the constrain t is still satisfied when u i = 0 due to the very large constant M [17]. (a) non-con ve x region (b) four con vex regions Fig. 2. Conv erting the non-con ve x region into con vex regions with geomet- rical relaxat ion Howe ver, whe n we want to d eploy the second U A V , an ad- ditional co nstraint ensurin g no overlappin g between cov erage areas, an d ther efore no ICI, is need ed. T o en su re this, the distance between the two U A Vs in the horizo ntal dim e nsion should be larger than 2 R . Theref ore, the placemen t of the second UA V is formulated as maximize x c 2 ,y c 2 ,u i X i ∈M u i (15) sub ject to ( x i − x c 2 ) 2 + ( y i − y c 2 ) 2 ≤ R 2 + M (1 − u i ) , ∀ i ∈ M ( x c 2 − x c 1 ) 2 + ( y c 2 − y c 1 ) 2 ≥ 4 R 2 u i ∈ { 0 , 1 } , ∀ i ∈ M where ( x c 2 , y c 2 ) is the horizon tal location of the secon d U A V . Unfortu n ately , the addition al con stra in t is non-convex which makes (15) extremely hard to solve. Altho ugh the integer variables in (13) can be addressed with advanced mixed integer progr amming techniqu es, u sing solvers suc h as MOSEK [17], the optim ization prob lem ( 15) which is a MI N L P prob lem with non-co n vex constra in t c a n not be straightfo rwardly solved. Even if we ap ply semidefinite relaxation ( SDR) tech niques to conv ert the quadr atic programs into the for m of semidefinite matrix which makes the no n-conve x constraint of (1 5 ) con vex, a problem with both po siti ve semid efinite matrix and integer variables is still un solvable with existing techniques [28]. In Fig. 2(a), the circle in white with radius R r epresents the coverage a r ea of the first aerial BS, an d the circle with radius 2 R represents th e area where there cannot be any placement of additional U A Vs witho ut inflicting ICI. Ac c ording ly , the region outside the gr een circ le with radiu s 2 R is the g eometrical representatio n of the non -conv ex con straint in (15). W e n otice that such a non -conve x region which spec ifies all the feasible locations o f the second U A V in the hor izontal dimension c an be di vided in to four linear region s which are co nvex. This is done b y ap p roximatin g th e gre e n circular area by a squ are area such tha t the orig inal gre en circle is surro unded by the square as illustrated in Fig. 2(b). Th erefore , instead of solving (15), we can solve fou r MINL P prob lems with different linear constraints, and each of the fou r prob lems has the following form maximize x c 2 ,y c 2 ,u i X i ∈M u i (16) sub ject to 5 Fig. 3. An example of feasible region definition , with two deplo yed aerial BSs, for the positio ning of the third BS ( x i − x c 2 ) 2 + ( y i − y c 2 ) 2 ≤ R 2 + M (1 − u i ) , ∀ i ∈ M y c 2 ≥ y c 1 + 2 R, if ( x c 2 , y c 2 ) ∈ A 1 x c 2 ≤ x c 1 − 2 R, if ( x c 2 , y c 2 ) ∈ A 2 y c 2 ≤ y c 1 − 2 R, if ( x c 2 , y c 2 ) ∈ A 3 x c 2 ≥ x c 1 + 2 R, if ( x c 2 , y c 2 ) ∈ A 4 u i ∈ { 0 , 1 } , ∀ i ∈ M where A 1 , A 2 , A 3 and A 4 are the four conv ex regions shown in Fig. 2 (b). The m aximum number of covered user s as well as the location o f the second U A V is then fo und amo ng the re su lts of f our MINL P pro b lems. Note that th e overlapping in the four conv ex regions will not affect th e final result an d is thus allowed. I f the optimal location of the second U A V is inside the overlapping a r ea, two of the fo ur optimization prob lems will give th e same solution which contains a maximum numb er of covered users. Howe ver, the effecti ve ar ea for placin g the second UA V is sligh tly reduced as a resu lt of co nsidering the infeasible r egion a s a square region. The op timization problem of placing th e k -th UA V ( k > 1 ) is for mulated as maximize x ck ,y ck ,u i X i ∈M u i (17) sub ject to ( x i − x ck ) 2 + ( y i − y ck ) 2 ≤ R 2 + M (1 − u i ) , ∀ i ∈ M ( x ck − x cj ) 2 + ( y ck − y cj ) 2 ≥ 4 R 2 , j = 1 , 2 , ..., k − 1 Algorithm 1 Algorithm for p lacing the k -th UA V with ge o - metrical relax ation Inputs: user lo cations, ( x i , y i ) ∈ M ; radius of coverage area, R ; locations of all dep loyed UA Vs ( x cj , y cj ) , j = 1 , 2 , ..., k − 1 Output: n umber of users covered by the k -th U A V , U k ; th e location o f the k -th U A V , ( x ck , y ck ) Initialization: j =1; z =1; m =0 . 1: while j < k do 2: conv erting the constraint ( x ck − x cj ) 2 + ( y ck − y cj ) 2 ≥ 4 R 2 into fo ur linear con straints which are x ck ≥ x cj + 2 R , x ck ≤ x cj − 2 R , y ck ≥ y cj + 2 R , an d y ck ≤ y cj − 2 R r espectively . 3: j = j + 1 . 4: end while 5: For each of the k − 1 UA Vs, one o f the four linear constraints is selected to form the intersection of these k − 1 regions. A total of 4 k − 1 intersections a r e gene rated and d e noted as C z , z = 1 , 2 , ..., 4 k − 1 . 6: while z < 4 k − 1 do 7: if C z = ∅ then 8: eliminate C z 9: else if C z ⊆ C q , q = 1 , 2 , ..., 4 k − 1 , q 6 = z then 10: eliminate C z 11: else 12: m = m + 1 . 13: C m k = C z . 14: obtain ( x m ck , y m ck ) , an d N m by solv in g (18) 15: end if 16: z = z + 1 . 17: e nd while 18: N max = max( N m ) , m = 1 , 2 ,...,N k M . 19: U k = N max . 20: ( x ck , y ck ) = ( x m ck , y m ck ) | N m = N max . u i ∈ { 0 , 1 } , ∀ i ∈ M where ( x ck , y ck ) and ( x cj , y cj ) den ote the ho rizontal locatio n of the k -th UA V an d the j -th U A V respe ctiv ely . F or each of the k − 1 non-conv ex constrain ts, g eometrical relaxatio n is utilized to co n vert it into four linea r co nstraints as illustrated above. As each linear constraint only specifies one of the four feasible region s gen erated b y a ce rtain UA V , the tr ue feasible region for the k -th U A V is the intersection of ( k − 1) sets, where each set is the union o f f our feasible region s. Evidently , the total feasible region which is n on-co nvex consists of se veral co nvex regions of rectan gular shape as sho wn in Fig. 3. The total n umber of c on vex regions depe nds on specific deployment but can be fo und throug h an elimin a tion meth od. T o be specific, for each of the k − 1 U A Vs, one of the four generated feasible r egions is selected and logic fun c tion is used to find th e intersection of these k − 1 selec ted regions to make sure the coverage area of the next aerial BS does n ot interfere with any pr eviously deployed aerial BSs. A total of 4 k − 1 intersections should be g enerated and we de n ote e a ch intersection as C z , z = 1 , 2 , ..., 4 k − 1 . After ob taining all the 4 k − 1 intersections, we eliminate all sets which ar e null sets, 6 i.e. C z = ∅ or sets which tu rns out to be sub sets of oth er generated sets, i.e. C z ⊆ C q , q = 1 , 2 , ..., 4 k − 1 , q 6 = z , an d the rem a in ing intersections are th e feasible regions we shou ld search for . An examp le is shown in Fig . 3. Her e, 5 feasible regions ar e found among 16 intersections for deploying the third aer ial BS after elimina ting all null sets and subsets. W e denote the total number of feasible regions f or deploying the k -th U A V as N k M . In this case, (17) can be r eformu lated as N k M MINLP problems, each h as the follo wing form maximize x m ck ,y m ck ,u i X i ∈M u i (18) sub ject to ( x i − x m ck ) 2 + ( y i − y m ck ) 2 ≤ R 2 + M (1 − u i ) , ∀ i ∈ M ( x m ck , y m ck ) ∈ C m k u i ∈ { 0 , 1 } , ∀ i ∈ M where C m k is the m -th feasible region of the k -th aerial BS, ( x m ck , y m ck ) is th e optimal location of th e k - th U A V in region C m k , m = 1 , 2 , ..., N k M . If we denote the numb er o f covered user s b y solving the m - th optimizatio n pr oblem as N m , an d den ote the maximum N m for all m a s N max , we hav e ( x ck , y ck ) = ( x m ck , y m ck ) | N m = N max . In addition, U k = N max , where U k denotes the numb er of covered users by the k - th aerial BS. For clarity , th e proposed geometrical relaxation method is summarized in Algorithm 1. I V . P RO P O S E D S I M U LT A N E O U S D E P L O Y M E N T M E T H O D W I T H K - M E A N S C L U S T E R I N G ( S D - K M ) The sho rtcoming o f the above prop osed algorith m is that it introdu c es exponentially increasing compu tational complexity due to the need to solve 4 k − 1 logic com bination operation s as well as multiple MINLP problems for finding the o ptimal location of th e k -th UA V , wh ich makes its use prohibitively complex when a large num ber of a erial BSs are need e d. As a result, ther e is a strong m otiv ation for a m echanism in wh ich multiple aeria l BSs ca n be deployed simu ltaneously without introdu c ing non- conv ex co nstraints. In this section, a method which simultaneously de p loys multiple U A Vs is proposed with the h elp of K -means clustering . K -me a n s clu stering, which is a well-known p artitional clustering method h a s bee n utilized in a variety of d isciplines [29]. In our particular scenario, we notice that the whole ta rget coverage area can b e divided into K subareas with b ounda r ies forming the V o ronoi diagram, by assigning user po in ts in to K clusters. The intellig e nt division of the target area brings great benefit to the d eployment of multiple U A Vs in sev eral senses. First, each subarea which is bou nded by f ew straig ht lines or line segments is a c o n vex region. W ithin each conve x region, we can solve an op timization problem similarly to (13) to find the best lo cation of a U A V so a max imum nu mber of user points within that subarea is enclosed. In addition , the bou ndary lines o f V oro noi diag ram en sure that the circles placed in each subarea will not overlap with each other, so the ICI is intrinsically av oided. Furthe rmore, a nu m ber of K U A Vs can be d e ployed simultaneo usly in correspo nding subareas, so the latency and d ependen ce o n pr eviously deployed aerial BSs with successi ve circle placement is solved. Last but no t the least, K -mean s clustering find s po tential clusterin g pr operties among user poin ts. The clustering p roperties in dicate a rea- sonable n umber of UA Vs to be deployed, hence a voiding the use of excessive U A Vs a nd saving bo th costs and power . The details of the prop osed me th od are illustrated in the following two subsection s. A. Applying K -mean s clustering and partitionin g the tar get ar e a W e assume the user set M contains a total of U tot users and we deno te arrays storing the location of user points by u i , wh ere u i = [ x i , y i ] , i = 1 , 2 , ..., U tot . T h e aim of apply in g K -me a n s clustering is o rganizing the U tot user points in to K clusters C k , k ∈ [1 , K ] , with the k -th cluster, k = 1 , 2 , ..., K , containing N k user points, ou t of which U k ≤ N k user po ints are covered. The par tition is very much based on a sum -of- squared- e rror criterio n [29] which is defin ed as e = K X k =1 U tot X i =1 u ki k u i − m k k 2 (19) where k . k denotes th e fr obenius nor m of a vector . u ki here is a Boolea n variable, indicating the state of i - th user of the k -th cluster . W e have u ki = 1 when ( x i , y i ) ∈ C k and u ki = 0 oth erwise. In addition , m k is an arra y storing th e center location of th e k -th cluster . Note that the center of a certain cluster is ob tained by calculatin g the mean value o f all user points classified into that cluster , which can be written as m k = [ m kx , m ky ] = ( 1 N k U tot X i =1 u ki x i , 1 N k U tot X i =1 u ki y i ) (20) where k = 1 , 2 , ..., K . Then the pro cedure of ap plying th e K -me a n s clu stering is concluded in four steps. 1) Random ly choo se K po ints in the target area as the center p oints and store them in m k , k = 1 , 2 , ..., K 2) Allocate eac h user point in M to the cluster with the closest cen ter C j , i.e., ( x i , y i ) ∈ C j when th e Euclidean distance between u i and the center of c lu ster j is smaller than the Euclid ean distan ce between u i and any other cluster center s. k u i − m j k < k u i − m k k , (21) i = 1 , 2 , ..., U tot , k = 1 , 2 , ..., K, j 6 = k 3) Recalculate the cluster cen ters as the m ean position o f all user points in each cluster . 4) Repeat th e above two steps u ntil there is no chang e for m k . Note that K -means clusterin g assumes the number o f clus- ters to g e nerate is k nown, which is not true fo r our case. Either excessi ve o r inad equate nu mber of gener ated subareas can affect the number o f users covered. The num ber of required clusters hig hly depen d s on user distributions, so variable K might be u tilized for various scen arios. For a ll cases, we first start with a max imum K value, denoted as K max , making sure adequate number of subareas ar e generated even for the case showing th e least clusterin g pr operty (unifo r m distribution). Then we try to r educe the number of K since user poin ts might 7 gather at some locations. As excessive pa r tition can split a single cluster into several parts, wh ich severely deterior ate the perfor mance o f the proposed method, we set a threshold d th indicating the minimu m allo wed distance between two cluster centers. If min( k m j − m k k ) < d th , j 6 = k , th is signifies that some o f the g enerated clusters are to o small as a re su lt of using too large K value. Then we re duce the value of K by one and reapply the K -means clustering. The above pr ocedure continues un til we h av e min( k m j − m k k ) ≥ d th . B. Solving th e o p timization pr o blem within each r e g ion After p artitioning the user poin ts into K cluster s an d , subsequen tly , di viding the who le target area into K subareas, we first n eed to find the largest allowed coverage area within each subarea to av o id ICI. As th e shape o f each sub area is a po lygon, it is likely th at c e rtain subareas can o nly accommo date circles with radii smaller than R . Assum e the k -th subarea is for med with S k line segments or straight lines, each line is expressed in the form of y = a kl x + b kl , l = 1 , 2 , ..., S k , where a kl and b kl represents the slope a n d offset of th e l -th b ound ary line of the k -th subarea respectively . It is known th a t for any poin t ( x d , y d ) , if y d − a kl x d − b kl < 0 , the po int is in the region b e low the line. On the contrary , if y d − a kl x d − b kl > 0 , the poin t is in the region ab ove the line. W e note that the boun d ary lines of each su barea also implicates a feasible region for p la c ing the circle c e nter as well as restricting the length o f radius. T o be specific, the distance b etween th e circle center and each b ounda ry line should be no smaller than th e length o f rad iu s of the c ircle. Therefo re, the r egio n for placing the circle can be fou nd by shifting the boundar y lines of each subarea. I f the clu ster center m k is in the region b elow a certain bound a ry line of the k -th subar e a, the correspo nding ne w line specifying the region for placing the c ir cle center can b e obtained by sh ifting the line downward along the y-ax is by L kl . Similarly , shifting the original boun dary line upward along the y-axis by L kl leads to the correspo nding new line when m k is in the region above the orig inal boun dary line. Here , L kl denotes the length to be shifted along the y-axis of the l -th bound ary line of the k -th subarea, and is calcu late d th rough L kl = R k max cos( | a kl | ) , k = 1 , 2 , ..., K, l = 1 , 2 , ..., S k (22) where R k max denotes th e maximum allowed rad ius o f th e circle placed in the k -th subar ea and | . | calculates the am plitude of a numbe r . Therefo re, R k max can b e found b y solv ing the following o ptimization pr o blem. maximize x ck ,y ck ,R k max R k max (23) sub ject to y ck − a kl x ck − b kl + L kl ≤ 0 , if m ky − a kl m kx − b kl ≤ 0 y ck − a kl x ck − b kl − L kl ≥ 0 , if m ky − a kl m kx − b kl ≥ 0 k = 1 , 2 , ..., K , l = 1 , 2 , ..., S k (a) (b) Fig. 4. The case for optimizing the radius in K -m eans circl e placement algorit hm: (a) flexibil ity in reaching addi tional users, (b) reducing po wer for a giv en user co ver age area . Algorithm 2 Iterative alg orithm fo r placing the k -th U A V Inputs: Initial radius R k ; an intermed iate v alue storing the change of radius, r it . Output: Set contain ing all covered user points, M k cov ; the location o f the k -th U A V , ( x ck , y ck ) ; th e op timal ra dius r k . Initialization: r it = 0 , r k = R k 1: while r it 6 = r k do 2: r it = r k 3: obtain ( x ck , y ck ) and M k cov by solv ing (2 4) and replac- ing R k with r it . 4: obtain r k by solv in g (25 ). 5: end while The radiu s of th e k -th circle is then R k = min( R, R k max ) . W ith fixed radii, we can find the optimal placem ent of cir c les for covering maximu m numbe r of user po ints within their correspo n ding suba r eas by solving the following p roblem maximize x ck ,y ck ,u i X i ∈M u i (24) sub ject to ( x i − x ck ) 2 + ( y i − y ck ) 2 ≤ R k 2 + M (1 − u i ) , ∀ i ∈ M y ck − a kl x ck − b kl + R k cos( | a kl | ) ≤ 0 , if m ky − a kl m kx − b kl ≤ 0 y ck − a kl x ck − b kl − R k cos( | a kl | ) ≥ 0 , if m ky − a kl m kx − b kl ≥ 0 u i ∈ { 0 , 1 } , ∀ i ∈ M k = 1 , 2 , ..., K, l = 1 , 2 , ..., S k The ab ove optimization p r oblem is a MINLP pr oblem without non-co n vex constraints and is near ly as easy to solve as (1 3). Note that with the help of K -means clustering, we on ly need to solve K optimization prob lems of this ty pe and K < K max is obtained in m o st cases as illustrated in th e first subsection. V . E N E R G Y E FFI C I E N T S I M U L TA N E O U S D E P L O Y M E N T M E T H O D W I T H V A R I A B L E R A D I U S ( S D - K M V R ) In the preceding section, a simultaneou s dep loyment method has been prop osed to take advantage of dividing the wh ole target area in to K co n vex sub areas. Howe ver, the prop osed 8 x-dimension (m) 0 500 1000 1500 2000 2500 3000 3500 y-dimension (m) 0 500 1000 1500 2000 2500 3000 3500 Fig. 5. SD-KM with inaccurate ULI method ba sed on K -mean s clustering can be further improved in terms of both the cov erage p robability and po wer ef ficie n cy . First, the ma x imum circle area do es n ot a lways cover a max- imum n umber o f user points in a irregu la r poly gon region , so variable rad ius m ay introdu c e furth er improved perfo rmance gain. As can be seen in Fig. 4(a) , user points ma y gath e r in a relativ e narr ow r egion where circles with large rad ii can not reach. More user points can th us be enclosed when a smaller coverage area is p laced. Furth ermore, the radii of coverage areas and hen ce th e transmit p ower of aerial BSs can be further reduced because there migh t be no user poin ts right on the bo rder of the circles. In Fig. 4 (b), origin al circle in red ob tained by the SK-KM in the previous section, can be shrunk into the circle in green which covers the same set of user po ints with a reduced tr ansmit power . In ord er to address both pr oblems at the same time, an iterativ e algor ith m is propo sed. W e assume aerial BS has a m inimum allowed coverage area with rad ius R min , th en the variable radius of th e k -th circle r k has a range of R min ≤ r k ≤ R k . W e first obta in the circle center ( x ck , y ck ) as well as the set con ta in ing th e covered u ser points, deno ted as M k cov with size U k by solv in g (24) with radius R k . Th en, we fix both ( x ck , y ck ) and M k cov , aim ing to find the minimum r k which can cover the same set of user p oints by solving the problem , minimize r k (25) sub ject to r k 2 ≥ ( x i − x ck ) 2 + ( y i − y ck ) 2 , ∀ i ∈ M k cov R min ≤ r k ≤ R k After solving r k , we r e place R k with r k and solve (24) again to find the update d user points enclo sed by the new circle. The above p rocedur e repeats until the radiu s r k does not ch ange anymore. Algor ithm 2 summaries the prop osed iterativ e alg orithm. In ad dition, as th e radii of th e K coverage Algorithm 3 Robust dep loyment of aerial BSs Inputs: Placement details ob tained from SD-KM o r SD- KMVR tech nique: radius o f coverage areas, R k ; ho rizon- tal location o f ae rial BSs, ( x ck , y ck ) ; location of covered user po ints, ( x i , y i ) , ∀ i ∈ M k cov Output: New horizon tal locatio n of aerial BSs, ( x ∗ ck , y ∗ ck ) ; new rad ius of c overage ar eas, R ∗ k . 1: Find th e minimum distance between ( x ck , y ck ) and the bound ary lin e s of the c o rrespon ding subar ea. 2: Obtain ( x ∗ ck , y ∗ ck ) by solving (31 ). 3: Calculate th e min im um distance between ( x ∗ ck , y ∗ ck ) a n d the bou ndary lines of the corresponding subarea. 4: Obtain R ∗ k from (32 ). areas a r e r educed and all the a erial BSs are alw ays placed at a position with fixed elev ation angle θ opt , the altitud e o f the k -th UA V can be fou nd by h k = r k tan( θ opt ) (26) Therefo re, the 3-D location of th e k -th a erial BS is marked as ( x ck , y ck , h k ) . Furtherm ore, the r e duced radii also reduce the path loss b etween U A Vs and users accordin g to ( 8) and thus r educe the required transmit power o f aerial BSs since we h av e P k t = P min + PL( r k ) (27) where P k t is the required tran smit power o f the k -th aerial BS and P min is the thr e shold receive power , below which the commun ication link is failed. The total required p ower o f the system is thus the sum o f transmit power of a ll aerial BSs which is a func tio n of bo th K an d r k , P total = K X k =1 P k t = K P min + K X k =1 PL( r k ) (28) Clearly , th e iterative alg orithm is more power efficient th an the prop o sed SD-KM algorithm at the cost of in creasing the computatio nal com plexity . V I . I M P E R F E C T U L I A N D R O B U S T D E P L OY M E N T In the pr a ctical U A V deployment the ULI m a y contain errors. Conseq uently , the coverage probability of the propo sed technique s may decrease drastically as a result of e stimating user locations inaccur ately . In this section, we for mulate a robust technique which is applicab le to bo th SD-KM and SD-KMVR to preserve the best coverage p erforma nce in the existence of inaccu rate ULI. W e model the estimated location of user i a s ( ∼ x i , ∼ y i ) = ( x i + e xi , y i + e y i ) , w h ere e xi and e y i are estimation error s following Gaussian distribution with zero mean an d standard deviation σ in meter s. For illustration, Fig. 5 shows an example d eployment of SD-KMVR in the existence o f imperfect ULI with σ = 50 m and L s = 3500 m. The small circles in red represen t real user location s, and SD-KMVR is ap plied based on estimated user locatio ns represented by do ts. It can b e seen that, user poin ts which are closer to th e centers of aerial BSs have better immunity to estimation errors. The coverage p r obability decreases when 9 User Density (users/km 2 ) 5 10 15 20 25 30 Average Execution Time (s) 0 0.05 0.1 0.15 0.2 (a) Number of Required Iterations 0 10 20 30 40 50 Cumulative Density Function 0 0.2 0.4 0.6 0.8 1 K-means SD-KMVR (b) Fig. 6. Computati onal c omplex ity: (a) avera ge exe cutio n ti me of solving a single MINLP problem by MOSEK solver , K = 1 ; (b) CDF of number of iterat ions requir ed for K -means clustering and SD-KMVR, K = 9 , λ s = 10 users / km 2 the users wh ich are consider ed to be covered ar e actually o ut of the coverage rang e o f th e corre sp onding a e rial BSs. It is clear that in creased robustness aga in st imperf ect ULI can be achieved by increasing the radii of coverage areas. If we assum e the maximum deviation between estimated location and real locatio n for any user po int is d th , where d th ≈ 3 σ , the performan ce lo ss of coverage pr o bability for the k -th aerial BS c a n be completely c ompensated wh en L k min = | R k − r ik | ≥ d th , ∀ i ∈ M k cov (29) where L k min denotes the minimum difference b etween r ik and R k of the k -th subarea. Th ere is clearly a trade-o ff between robustness and required tran smit power with regard to the radius. Theref ore, the aim of the robust design, wh ich is maximizing the number o f covered user points in th e existence of inaccurate ULI, is equi valent to max imizing L k min with m in- imum tr a nsmit p ower . W e no te that the u ser p oints are usually distributed unevenly within the c orrespon ding coverage area for both SD-KM and SD-KMVR techniqu es. This causes a part of user po in ts having a much larger gro und distance to the aerial BS than the r est of u ser points. Therefore, shiftin g the circle center to min imize the ma x imum grou nd d istance between user points covered b y the k - th aerial BS and the k - th circle cen te r can red u ce the resulting rad ius and th us redu ce the required tran sm it power . W e d enote th e d istance between ( x ck , y ck ) and th e boundary lin es of th e co rrespond ing subarea as d kl , l = 1 , 2 , ..., S k . Then th e minimu m v alue among all the S k distances is d k min = min( d kl ) . In order to a void ICI, the horizon tal center of th e k - th ae r ial BS can only move within a circular ar ea with radiu s d k min . Ther e fore, the correspo nding optimization p roblem is formulated a s minimize x ∗ ck ,y ∗ ck max i ∈ M k cov ( q ( x ∗ ck − x i ) 2 + ( y ∗ ck − y i ) 2 ) (3 0) sub ject to q ( x ∗ ck − x ck ) 2 + ( y ∗ ck − y ck ) 2 ≤ d k min k = 1 , 2 , ..., K Here, ( x ∗ ck , y ∗ ck ) is the new ho rizontal ce nter of the k - th aerial BS to optimize, based on the center co ordinates ( x ck , y ck ) ob- tained by either SD-KM o r SD- KMVR. Note that the o bjective function of the above o p timization problem im plicitly includes the constraint that all the o riginally covered user points are st ill covered. (30) is eq uiv alen t to minimizing an auxiliary variable d k representin g the maximu m grou n d distance between k - th aerial BS and covered u ser poin ts accor ding to minimize x ∗ ck ,y ∗ ck d k (31) sub ject to q ( x ∗ ck − x i ) 2 + ( y ∗ ck − y i ) 2 ≤ d k , i ∈ M k cov q ( x ∗ ck − x ck ) 2 + ( y ∗ ck − y ck ) 2 ≤ d k min k = 1 , 2 , ..., K After obtaining the new cen ter location ( x ∗ ck , y ∗ ck ) , we recalcu- late the m inimum gr ound distance between the k - th aerial BS and the correspon ding b ound a ry lines and denote it as d k ∗ min . The max imum allowed r a dius within the k -th sub area is then R k ∗ max = min( R, R k + d k ∗ min ) . Theref ore, th e radiu s of the k -th coverage area R ∗ k is R ∗ k =  d k + d th , if d k + d th ≤ R k ∗ max R k ∗ max , if d k + d th > R k ∗ max (32) Note tha t the rad iu s R ∗ k is no t necessarily larger than R k , especially when the technique is applied to SD-KM. Therefor e, a redu ced transmit power is sometimes obtaine d. F or clarity , the above proced ure is summarized in Algorithm 3. V I I . C O M P U TA T I O NA L C O M P L E X I T Y A NA L Y S I S In this sectio n, we stud y the computatio n al complexity of the proposed techniques in term s of flo a tin g-poin t opera tions required . Following [ 30], [3 1], th e compu tational costs are calculated based on real-valued addition s, sub tractions, multi- plications, divisions and comparisons. A. Complexity o f SD- GR For deploying the k - th aerial BS ( k > 2 ), we first need to find a total of 4 k − 1 candidate sets. Eac h of the 4 k − 1 sets is an in tersection of k − 1 sets, and fo rming e a ch in tersection needs 4( k − 2) comparisons in the worst case for both x and y dimensions. Therefore, the complexity of finding th e candidate regions for all the K aer ial BSs needs K P k =3 2( k − 2 )4 k floating- point o perations, which can b e simp lified as C 1 GR = O { 16 K − 14 9 · 4 K +1 } (33) In the eliminatio n pr o cess, the complexity arises fro m the search fo r feasible regions from all can didate sets. There are a total of K P k =2 4 K − 1 candidate sets, and the resulting complexity is ob ta in ed as C 2 GR = O { 1 3 · 4 K +1 } (3 4) The a verage complexity of solv in g (18) does n ot have a closed form solution, an d it is in general difficult to determine . Sin c e this co mplexity is in volved in all the dep loyment schemes, we denote this a s C MINLP , and represent the comp lexity of all 10 T ABLE I C O M P U TATI O N A L C O M P L E X I T Y O F T H E P RO P O S E D T E C H N I Q U E S Method Computati onal costs SD-GR O{ 6 K − 11 9 · 4 K +1 } + K C MINLP SD-KM O { N KM K U tot n it + K P k =1 (3 + S k ) 1 . 5 } + K C MINLP SD-KMVR O { N KM K U tot n it + K P k =1 (3 + S k ) 1 . 5 + n vr it K P k =1 ( U k + 1) 1 . 5 } + ( n vr it + 1) K C MINLP Robust O { N KM K U tot n it + K P k =1 (3 + S k ) 1 . 5 } SD-KM + K C MINLP + O { K P k =1 U k + K 2 } Robust O { N KM K U tot n it + K P k =1 (3 + S k ) 1 . 5 SD-KMVR + n vr it K P k =1 ( U k + 1) 1 . 5 } + ( n vr it + 1) K C MINLP + O { K P k =1 U k + K 2 } schemes as a function of th is com plexity . Therefo re, the total computatio nal com plexity for SD-GR technique is C GR = C 1 GR + C 2 GR + K C MINLP (35) = O{ 6 K − 1 1 9 · 4 K +1 } + K C MINLP T o characte r ize the co mplexity of solving a single MINLP , instead o f an analy tical expression, we emp loy the average execution time against various user density as shown in Fig. 6 ( a). B. Complexity of S D - KM W e the n consider the co mplexity of K - means clustering . Here, we denote th e average numbe r of iterations until con- vergence as n it . T he cumulative distribution functio n (CDF) describing the conv ergence for K -means clustering is shown as the d ashed lin e in Fig. 6 (b ). Within each iteration, the prop osed sch e me inv olves th ree steps. The first stage calculates the Euclid ean distance between each user po int and cluster cen ters, wh ich in volves two subtr actions, two multiplications, o ne additio n and one square r oot. For U tot user points and K clusters, calculating all Euclidean distances requires O{ 6 K U tot } floating-po int operations. The secon d stage allocates each u ser point to the cluster with the closest center, wh ich takes O { U tot ( K − 1) } c o mparison s. Further- more, we need to recalcu la te the cluster centers fo llowing ( 2 0), which includes U tot multiplications, U tot − 1 additions and one division for bo th x an d y d imensions. The n the costs for the third stage is O { 4 K U tot } . Th erefore, the resu lting computatio nal com plexity of k - means clustering is O{ n it (6 K U tot + U tot ( K − 1 ) + 4 K U tot ) } ≈ O { K U tot n it } (36) 0 1000 2000 3000 y-dimension (m) 0 500 1000 1500 2000 2500 3000 CPT 0 1000 2000 3000 0 500 1000 1500 2000 2500 3000 SD-GR x-dimension (m) 0 1000 2000 3000 y-dimension (m) 0 500 1000 1500 2000 2500 3000 SD-KM x-dimension (m) 0 1000 2000 3000 0 500 1000 1500 2000 2500 3000 SD-KMVR Fig. 7. Aerial BS placement with proposed tec hniques T o o btain a r e a sonable value of K , the pr oposed a l- gorithm repeats the ab ove steps for N KM times until min( k m j − m k k ) ≥ d th as shown in Section IV , an d we have C 1 KM = O{ N KM K U tot n it } (37) Moreover , we need to calculate the maximum allowed r adius within each cluster according to (23). F ollowing [32], (23) is a conve x problem solved by interior-point m ethods with the following co mplexity: C IP = O{ ( E + F ) 1 . 5 E 2 } (38) where E is the n umber o f variables, and F is the numb er of constraints in an optimiza tion pr oblem. As we h av e a total of S k constraints and 3 variables for solving th e k - th optimization pr oblem, the comp utational costs o f finding all maximum allowed radius is C 2 KM ≈ O { K X k =1 (3 + S k ) 1 . 5 } (39) The computational c o mplexity of solving (24) is again approx- imated b y C MINLP . Ac c o rdingly , the total com putationa l cost of SD-KM is C KM = C 1 KM + C 2 KM + K C MINLP (40) = O{ N KM K U tot n it + K X k =1 (3 + S k ) 1 . 5 } + K C MINLP C. Comp lexity o f S D-KMVR The SD-KMVR technique is based on SD-KM techniq ue and thus in volves all operation costs of SD-KM. I n addition, SD-KMVR is an iterative algorithm, an d deno te the a verage number of iter a tio ns requ ired as n vr it . CDF of the numbe r of 11 (a) HPP IPP PCP Coverage probability 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 CPT SD-GR SD-KM SD-KMVR (b) HPP IPP PCP 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 CPT-error SD-GR-error SD-KM-error SD-KMVR-error (c) HPP IPP PCP 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 SD-KM-error Robust SD-KM SD-KMVR-error Robust SD-KMVR Fig. 8. User-c ov erage probability for dif ferent types of user distrib ution: (a) with perfect UL I, (b) with imperfec t UL I, (c) with rob ust technique , K =4 required iter ations for SD-KMVR tech nique is also shown in Fig. 6 (b ). W ithin each iteration, the costs of solv ing (24) is C MINLP and the c o sts of solving ( 2 5) is O { K P k =1 ( U k + 1) 1 . 5 } , since the number of c o nstraints of ( 25) is U k . Therefore, the overall costs of SD-KMVR tec h nique is C KMVR = C KM + O{ n vr it K X k =1 ( U k + 1) 1 . 5 } (41) + n vr it K C MINLP = O{ N KM K U tot n it + K X k =1 (1 + S k ) 1 . 5 + n vr it K X k =1 ( U k + 1) 1 . 5 } + ( n vr it + 1) K C MINLP D. Complexity of Ro bust T e chn iq ue The additiona l com putational co sts of apply in g the ro bust technique co me fro m (30 ), which is again solved b y interio r- point metho d. Based on (30), E = 2 and F = U k + K − 1 , which lead s to C robust ≈ O{ K X k =1 U k + K 2 } (42 ) For clarity , the co mputation al complexity of th e prop o sed technique s is summarized in T ab le I. Note that the benchmar k CPT is origin ally designed for maximizin g the c overage ar ea instead of m a ximizing the nu mber of users covered, so the location o f aerial BSs is fixed fo r a specific target area wh en R is g iv en, and the tec hnique has n egligible com plexity . T o compleme n t the above complexity analysis, in Fig. 11 in the following we sho w a co mplexity co mparison of th e prop osed SD-KM, SD-KMVR and robust techniques in term s of average execution time. V I I I . S I M U L A T I O N R E S U LT S A N D A N A L Y S I S In this section, we assume multiple aerial BSs are vertically deployed at a position wh ic h m aximizes the coverage area in Number of UAVs 5 10 15 Coverage probability 0 0.2 0.4 0.6 0.8 1 HPP Number of UAVs 5 10 15 0 0.2 0.4 0.6 0.8 1 IPP SD-GR SD-KM SD-KMVR CPT Number of UAVs 5 10 15 0 0.2 0.4 0.6 0.8 1 PCP Fig. 9. User-co ver age probabi lity ve rsus number of U A Vs deplo yed, K =16, 15 and 10 for HPP , IPP and PCP correspondingl y T ABLE II S I M U L AT I O N PA R A M E T E R S paramete r value a 9.61 b 0.16 η LoS 1 η NLoS 20 f c 2.5Ghz c 3 · 10 8 m/s θ opt 42 . 44 ◦ an urban environmen t. Ther e f ore, with the parameters shown in T able II, th e rad ius R corr e sponds to a path loss threshold PL th = 100 dB is calculated as R = 7 07 m. The minimum radius o f each coverage area and the m in imum allowed distance between clu sters are assumed to be R min = R 2 and d th = R 2 respectively . The value of K max is set equal to the num ber of circles r e su lting from CPT . W e assume P min = − 7 0 dBm and all the a erial BSs initially have the same transmit power P t = 30 dBm. Th r ee different spatial point pro cesses ar e utilized for m odeling the user distribu- tion, which ar e Homo geneou s Po isson process ( HPP) with λ s = 5 user s / km 2 , Inh omog e neous Poisson p rocess (IPP), with λ ( x, y ) = 5( x 2 + y 2 ) use rs / km 2 and Poisson cluster process (PCP) re spectiv ely . Specifically , the ’p arent’ points of cluster p rocess are gen e r ated following HPP with λ p = 1 users / km 2 and the ’children’ po in ts are ge n erated with λ c ( x, y ) = α 2 π σ 2 e − 1 2 σ 2 ( x 2 + y 2 ) (43) where α = 0 . 9 and σ = 0 . 02 . T o e valuate the benefit of the propo sed a lg orithms, nu merical results based on Monte Carlo simulations o f th e propo sed SD-GR, SD-KM , SD-KMVR and the ro bust tech n iques are compared w ith the p erform ance of CPT which serves as the ben c hmark. Note that due to the exponentially increa sin g comp utational complexity o f SD-GR as shown in Section VI I, the maximum K value we use for the SD-GR algorith m is fo ur . The horizon tal center of all d e p loyed UA Vs m ust fall in sid e the target area, and we assum e th e coverage areas outside th e 12 Side length of target area (meter) 1000 2000 3000 4000 5000 Number of UAVs required 0 2 4 6 8 10 12 14 16 CPT SD-KM/SD-KMVR/Robust-HPP SD-KM/SD-KMVR/Robust-IPP SD-KM/SD-KMVR/Robust-PCP Side length of target area (meter) 1000 2000 3000 4000 5000 Total transmit power (dBm) 28 30 32 34 36 38 40 42 CPT SD-KM-HPP SD-KMVR-HPP Robust SD-KM-HPP Robust SD-KMVR-HPP SD-KM-IPP SD-KMVR-IPP Robust SD-KM-IPP Robust SD-KMVR-IPP SD-KM-PCP SD-KMVR-PCP Robust SD-KM-PCP Robust SD-KMVR-PCP Fig. 10. Required number of a erial BSs and total transmit po wer versus size of target area target area will no t cau se fu rther in terferenc e to u sers o u tside the in terested region. T o illustrate the function of the p roposed solutions, example drone placement distributions are shown ba sed on simulation for th e benchmark CPT , SD-GR, SD-KM, and SD-KMVR in Fig. 7, assuming a PCP distribution of users and L s = 3 km. The benchmar k CPT simp ly places circles with same size in a way that maximum coverag e is achiev ed an d no ne of these circles overlap. For CPT , the number of circles to be placed in a squ are target area N cp depend s on th e size of target area, which is N cp =  L s 2 R  2 where ⌈ . ⌉ is a ceiling functio n . A. Coverag e Pr oba bility Intuitively , it can be seen that, coverage probability highly depend s on th e user distrib ution when CPT is applied. In addi- tion, SD-GR method always aims to cover th e most number of users in the rema ining region b u t the placement of the U A Vs is restricted by the previously deployed BSs, which limits th e achiev able perfo rmance. Th e SD- KM method aims to find the clustering pro perties among user points an d is thu s more robust to the chang e of user distributions. Howe ver, the shap e of subareas limits the movement of U A Vs within each V oronoi cell, which may lead to users gatherin g in a relative narrow region of the subarea un covered. This drawback is solved by applying SD-KMVR, which also red uces the tran smit power to a large e xtent. T he above effects, are captured in Fig. 8 (a), which illustrates the achieved user-coverage pr obability , for d ifferent typ es of user distributions. For a fair compar ison of the achiev able coverage probab ility , a target ar ea with L s = 4 R is assumed , within wh ich all fo u r method s ca n horizon tally dep loy a maximum of fou r circles. It can be observed that th e coverage pro bability of all techniqu es are affected by user distributions. Note that th e perfo rmance of CPT decr eases while the achieved coverage prob ability of all the other techn iques incr ease wh e n th e u ser p oints te n d to have an u neven distribution, especially when clusters are form ed. Specifically , the propo sed SD-KMVR techniqu e a c hieves an up to 30% higher coverage proba b ility than th e be n chmark . The result is expected, b e cause CPT p lace circ le s in fixed locations for a giv en target area no matter how the user s are distributed, wh ic h highly deteriorate the cov erage perfo rmance when clusters are formed ou tside the coverage areas. On the contrary , our proposed algorithms ar e not restricted to fixed lo c ations, but in stead can be flexibly placed according to the cha nge of user distribution. When hetero geneity of user d istribution is introdu ced, especially when clusters are formed , user points are located clo ser to eac h o ther an d there is cor respond in gly a better c h ance to cover mor e users with in each ap plied circles. Mo re th a n 90 per centage of users ar e covered by applyin g the prop osed tech n iques when user s are distributed fo llowing PCP . The coverage perform a nce of the p roposed techn iq ues with imperfect ULI is shown in Figur e. 8 (b). It can b e seen that the performan ce of all techn iques except CPT d ecreases w h en introdu c ing imperf e ct ULI . T he per f ormanc e of CPT rem ains unchan ged b ecause the plac e m ent rule o f CPT is irr elev an t to ULI. Note tha t SD-KM metho d shows much better immunity to im perfect U L I th an SD-KMVR. This is as expected, since SD-KM m ethod utilizes a larger coverage area causing a larger distance be twe e n user points and the border of circles than SD- KMVR. Th e perfo r mance loss is greatly compensated when the propo sed robust technique is applied as sh own in Figure. 8 (c). The increased coverage probability is ac hieved as a resu lt of in creasing coverage area as well as relocating the aer ial BSs. In real scenario s, o nly a limited num ber of U A Vs may be av a ilable fo r deploymen t. Therefore, it is me a ningfu l to examine the coverage p robability of the propo sed tech niques versus the numb er of av ailable ae rial BSs. W e assume L s = 5 km, the K value used for HPP , IPP and PCP ar e 16, 15 and 10 respectively , wh ich are the a verag e K values we n eed for target areas of this size, and the corr e sp onding resu lts are shown in Fig. 9. As expected , n o matter what user d istribution is, the proposed SD-GR technique significantly o utperfo rms other technique s since the UA Vs ar e always d eployed in a position such that a maximu m n u mber of remainin g user points are covered. Moreover , it can be seen tha t SD-KMVR achieves an up to 10% pe r forman ce g ain compared to SD-KM. Th e SD- KM, SD-KMVR and the CPT tech n iques have compar a ble perfor mance w h en users ar e distributed uniformly . This is because when u ser s ar e un iformly distributed, th e K -mean s clustering metho d will di vide the target ar ea in a similar way as we using CPT . Whe n users are distributed fo llowing a non-u niform distribution, CPT o utperfo rms the proposed SD- KM and SD-KMVR meth ods when o nly a small num ber o f U A Vs are av ailable. In th is ca se, a similar num b er of U A Vs are deployed in a mor e tight way than the c ircle packing technique , causin g reduce d coverage ar eas of certain aerial BSs and hence a smaller nu mber of covered user s within th ese areas. When CPT is utilized, circles placed at position s wher e the user points are densely located ca n thus cover more user points due to the larger coverage area. Ho wev er , when clusters are pr esent, the propo sed SD-KM an d SD-KMVR techniq ue regain the superio rity as the U A Vs are d e p loyed at positions where clu stering prop erties ar e f ound. 13 User Density (users/km 2 ) 5 10 15 20 25 30 Average Execution Time (s) 0 5 10 15 20 25 30 SD-KM Robust SD-KM SD-KMVR Robust SD-KMVR (a) Number of Aerial BSs 5 7 9 11 13 15 Average Execution Time (s) 0 5 10 15 20 25 30 SD-KM Robust SD-KM SD-KMVR Robust SD-KMVR (b) Fig. 11. A verage exec ution time for the proposed techniqu es: (a) versus v arious user densit y , K =9 (b) versus v arious number of aerial BSs, λ s = 5 users / km 2 B. Energy Efficiency In Fig . 10 , we co mpare th e nu mber of requir ed aerial BSs and the req uired total tran smit power for SD-KM, SD- KMVR, Robust SD-KM, Rob u st SD-KMVR and CPT . It c an be seen that, SD-KM and SD-KMVR as well as their robust technique s require a smaller numbe r of UA Vs comp ared to CPT when u sers are n ot un iformly distributed in the ta rget area. Note that the gap between the pr oposed algo rithms and CPT becomes larger as the size of the target area increases, with SD-KM and SD-K M VR red ucing the n umber o f U A Vs required d own to 60% of that f or CPT . An up to 15% po wer is sa ved by applying SD-KMVR wh en users are distributed following PCP . In addition, tho ugh SD-KM and SD-KMVR require th e same number of aerial BSs to be deployed in order to max imize the coverage pro b ability , the SD- KMVR technique saves u p to 10 % power as can be seen. Robust SD- KMVR c o nsumes app roximately 1% mo re power than SD- KMVR as a result o f increasing the coverage area, but a c le a rly reduced transmit po wer is still obtained co mpared to SD-KM technique . It is worth highlighting that, Robust SD-KM is e ven more power efficient than SD-KM, wh ic h mean s the min imum distance between covered user points and the bo rder of the correspo n ding coverage area after reloc a ting the circle cen ter is larger than d th in most cases. I ndeed, the redu ced number of U A Vs saves operation costs and the r educed transmit power can pr o long the operation time of aerial BSs. C. Computa tional Comp lexity In Fig. 11, we characterize the comp lexity of SD-KM, SD- KMVR and their robust techniques in terms of the average execution time. The u ser p oints are distributed following HPP and an I n tel Cor e i7-67 00 2.6 G H z CPU computer is used for perfor ming th e simulation. The average ex ecution time versus various user densities is shown in Fig. 1 1(a) with K =9, while the average execution tim e versus various numb er of aer ial BSs is presented in Fig. 11 (b) with λ s = 5 users / km 2 . I t can be ob served from bo th figure th at, SD-KMVR technique takes more execution time than SD-KM, a n d the complexity of SD- KMVR in creases more faster . T o be specific, the execution time of SD-KM VR increases appr oximately 40% and 10 5% faster than SD-KM, as the num ber o f user poin ts an d the number of aerial BSs increase correspond ingly . Similar trends is found for Robust SD-KMVR a nd Robust SD-KM. Note that the ro bust techniqu es also introduc e increased compu ta tio nal complexity , which is consistent with the analytic results shown in Section VII. I X . C O N C L U S I O N In this paper, we have stud ied the efficient de p loyment of multiple aerial BSs in or der to maximize the number of covered users while avoiding ICI. A suc c essi ve dep loyment method co n verting each no n-conve x constra in t into four linear constraints is firstly prop osed with geom etrical relax ation. In order to red uce the computatio nal com plexity , we then propo se a simultaneo us d eployment metho d based o n K - means clustering, which con verts the target area into K convex subareas an d find the optimal location of U A Vs within each subarea. Mor eover , an iter a tive algorithm is utilized to re d uce the required transmit po wer while further improve the coverage probab ility . T o increase the robustness again st imperfec t ULI, a robust tech nique wh ich relocates the aerial BSs befor e increasing the radius of coverage areas is prop osed. Simulation results show th at th e prop osed algorith ms increase the cover- age pe r forman ce b y up to 30%. In ad dition, SD-GR m ethod is suitable fo r scenario s where a sma ll nu mber of U A Vs are av ailab le, and the iterativ e algo rithm achieves an up to 15% improved power efficiency at a co st of increased computatio n al complexity . Performan c e loss in the existence of inaccu rate ULI is also greatly compensated with the pro posed robust technique . R E F E R E N C E S [1] K. Namuduri, Y . W an, and M. Gomathisank aran, “Mobile Ad Hoc Networ ks in the Sky: State of the Art, Opportunities, and Challenges, ” in Pr oceedings of the Secon d ACM Mobi Hoc W orksho p on Airborne Network s and Communica tions , ser . ANC ’13. New Y ork, NY , USA: A CM, 2013, pp. 25–28. [2] E. W . Frew and T . X. Bro wn, “Airborne Communicati on Networks for Small Unmanned Aircraft Systems, ” Pr oceedi ngs of the IEEE , vol. 96, no. 12, pp. 2008–2027, Dec 2008. [3] A. Merw aday and I. Guven c, “U A V Assisted H eterogeneous Network s for Public Safety Communications, ” in 2015 IEEE W ire less Communi- cations and Netw orking Confer ence W orkshops (WCNCW) , March 2015, pp. 329–334. [4] Y . Zeng, R. Z hang, and T . J. Lim, “W ireless Communicati ons with Unmanned Aeria l V ehicles: Opportuniti es and Challeng es, ” IE EE Com- municati ons Magazi ne , vo l. 54, no. 5, pp. 36–42, May 2016. [5] C. M. Cheng et al. , “Maximizing Throughput of UA V -Rela ying Net- works with the Load-Carry-and-De li ver Paradi gm, ” in 200 7 IEEE W ire- less Communic ations and Network ing Confere nce , March 2007, pp. 4417–4424. [6] Y . Zeng , R. Zhang, and T . J. Lim, “Throughp ut Maxi mizati on for U A V - Enabled Mobile Rel aying Systems, ” IEEE T ransactions on Co mmuni- cations , vol. 64, no. 12, pp. 4983–4996, Dec 2016. [7] S. Zhang, H. Zhang, Q. He, K. Bian, and L. Song, “Joint Traj ecto ry a nd Po wer Optimizati on for U A V Relay Networks, ” IEE E Communicat ions Letter s , vol . 22, no. 1, pp. 161–164, J an 2018. [8] M. Hua, C. Li, Y . Huang, and L . Y ang, “Throughput m aximization for U A V -enabl ed wireless power transfer in relaying s ystem, ” in 2017 9th Internati onal Confer ence on WCSP , Oct 2017, pp. 1–5. [9] B. Pearre a nd T . X. Brown, “Model-free Traj ectory Opti mizati on for W irel ess Data Ferries among Multiple Sources, ” in 2010 IEEE Globecom W orksho ps , Dec 2010, pp. 1793–1798. 14 [10] Y . Zeng and R. Z hang, “Ener gy-Ef ficien t U A V Communicati on W ith Tra jecto ry Optimiza tion, ” IEEE T ransact ions on W irele ss Communica- tions , vol. 16, no. 6, pp. 3747–3760, June 2017. [11] E. Kalanta ri, I. Bor-Y aliniz, A. Y ongacog lu, and H. Y ani komero glu, “User Ass ociation and Bandwi dth Alloca tion for T errestr ial and Aerial Base Stations with Backhaul Considerati ons, ” in 2017 IE EE 28th An- nual International Symposium on P ersonal, Indoor , and Mobile Radio Communicat ions (PIMRC) , Oct 2017, pp. 1–6. [12] B. Galkin, J. Kibilda , and L. A. DaSilv a, “Deployment of U A V -mounted Access Poin ts According to Spatial User Loca tions in T w o-tier Cellular Networ ks, ” in 2016 W irele ss Days (WD) , March 2016, pp. 1–6. [13] J. Xu, Y . Zeng, and R. Zhang, “U A V -enabled Multiuser W irel ess Po wer Tra nsfer: Traj ectory Design and Energy Optimi zatio n, ” in 2017 23rd Asia-P acific Confere nce on Communications (APCC) , Dec 2017, pp. 1– 6. [14] D. Y ang, Q. Wu, Y . Zeng, and R. Z hang, “Energy T rade-of f in Groun d- to-U A V Communication via Traject ory Design, ” IEEE T ransaction s on V e hicula r T ec hnology , pp. 1–1, 2018. [15] J. Lu, S. W an, X. Chen, and P . Fan, “Energy-Ef ficient 3D U A V -BS Placemen t versus Mobile Users’ Density and Ci rcuit Powe r, ” in 2017 IEEE Globecom W orkshops (GC Wkshps) , Dec 2017, pp. 1–6. [16] L. Gupta, R. Jain, and G. V a szkun, “Survey of Im portant Issues in U A V Communicat ion Network s, ” IEE E Communicat ions Surve ys T utorials , vol. 18, no. 2, pp. 1123–1152 , Secondqua rter 2016. [17] M. Alzenad, A. El-Ke yi, F . Lagum, and H. Y ani komero glu, “3-D Placemen t of an Unmanned Aerial V ehi cle Base Station (UA V -BS) for Energy-Ef ficient Maximal Cov erage, ” IEE E W ireless Communi cation s Letter s , vol . 6, no. 4, pp. 434–437, Aug 2017. [18] A. Al-Houra ni, S. Kandee pan, and S. Lardner , “Optimal LAP Altitude for Maximum Cov erage, ” IEEE W irele ss Communications Letters , vol. 3, no. 6, pp. 569–572, Dec 2014. [19] M. Alzenad et al. , “3-D Placeme nt of an Unmanned Aerial V ehi cle Base Stat ion for Maximum Cov erage of Users Wi th Dif ferent QoS Require ments, ” IEEE W ireless Communicatio ns Letters , vol. 7, no. 1, pp. 38–41, Feb 2018. [20] R. I. Bor- Y aliniz, A. El-Ke yi, and H. Y anik omeroglu, “Efficie nt 3- D Placement of an Aerial Base Station in Next Generation Cellular Networ ks, ” in 2016 IEEE International Con fer ence on Communi cations (ICC) , May 2016, pp. 1–5. [21] E. Kalantari, M. Z. Shaki r , H. Y anik omeroglu, and A. Y ongac oglu, “Backh aul-a ware Robust 3D Drone Placeme nt in 5G+ W ireless Net- works, ” in 2017 IEEE Internati onal Confe re nce on Communications W orksho ps (ICC W orkshop s) , May 2017, pp. 109–114. [22] M. Mozaf fari , W . Saad, M. Bennis, and M. Debbah, “Drone Small Cell s in the Cl ouds: Design, Deployment a nd Performance Ana lysis, ” in 2015 IEEE Globa l Communic ations Confer ence (GLOBE COM) , Dec 2015, pp. 1–6. [23] ——, “ Effic ient Depl oyment of Multi ple Unmanned Aerial V ehicles for Optimal W irel ess Cov erage , ” IEEE Communicati ons Letter s , vol. 20, no. 8, pp. 1647–1650, Aug 2016. [24] J. L yu, Y . Zeng, R. Zhang, and T . J. Lim, “Placemen t Optimizati on of U A V -Mounted Mobile Ba se Stat ions, ” IEEE Communic ations Letters , vol. 21, no. 3, pp. 604–607, March 2017. [25] D. Xu and Y . Tian, “A Comprehensi ve Surve y of Clustering Algo- rithms, ” Annals of Data Science , vol. 2, no. 2, pp. 165–193, Jun 2015. [26] J. Illian, P . Pentti nen, H. Stoyan, a nd D. St oyan , Statistical Analy sis an d Modell ing of Spatial P oint P atterns , ser . Statistics in Practice . Wil ey , 2008. [27] A. Baddele y , Spatial P oint Pr ocesses and their Applicat ions . Nether- lands: Springer , 2007, vo l. 1892, pp. 1–75. [28] Z. q. Luo, W . k. Ma, A. M. c. So, Y . Y e , and S. Zhang, “Semide finite Re- laxat ion of Quadrati c Optimiza tion P roblems, ” IEEE Signal Proce ssing Magazi ne , vo l. 27, no. 3, pp. 20–34, May 2010. [29] R. Xu and D. Wunsc h, Clustering . W iley-IEEE Press, 2009. [30] S. Bo yd and L. V a ndenber ghe, Conv ex Optimization . New Y ork, NY , USA: Cambridge Unive rsity P ress, 2004. [31] F . Liu, C. Masouros, A. Li, H. Sun, and L. Hanzo, “MU-MIMO Commu- nicat ions W ith MIMO Radar: From Co-Existe nce to Joint Transmission, ” IEEE Tr ansactions on W irel ess Communicati ons , v ol. 17, no. 4, pp. 2755–2770, April 2018. [32] A. Li, C. Masouros, F . Liu, and A. L. S w indlehurst, “ Massi ve MIMO 1- Bit DA C T ransmission: A Lo w-Comple xity Symbol Sc aling Approach, ” 2017.