I Introduction
To support the exponentially growing data traffic and number of devices in future wireless network, nonorthogonal multiple access (NOMA) was proposed to serve multiple users simultaneously on a same resource block [1, 2, 3]. The key idea of NOMA is to combine the superposition coding at transmitter and successive interference cancellation (SIC) at receiver, as such the spectral efficiency of wireless systems can be significantly improved [4, 5]. Due to its capability to provide superior spectral efficiency and massive connectivity, NOMA has been applied to many aspects of wireless communication systems, e.g., the Internet of Things network [6], and the ultrareliable and lowlatency communications network [7]. User clustering, also named as user grouping/pairing, is one fundamental issue of NOMA. The impacts of user clustering on the system performance have been intensively investigated from both performance analysis and system design perspectives. For example, [8]
studied the influences of user pairing on the outage probabilities of users in the NOMA system.
[9]proposed to use a branchandbound based algorithm to solve the joint user pairing and power allocation optimally with the worst case computation complexity of NPhard. To develop a computationefficient algorithm for the user clustering problem, matching theory based heuristic algorithms were presented in
[9, 10, 11]. Recently, to further reduce the computational complexity, [9] proposed to use the kmeans based approach to perform unsupervised user clustering by exploring the users’ position information.Recently, unmanned aerial vehicles (UAVs) assisted wireless communications has received considerable attentions due to its advantages to provide real time and high throughput services [12]. Compared to the conventional terrestrial wireless communication systems, the development of UAV has also created a fundamental paradigm shift to facilitate fast and highly flexible deployment of communication infrastructures. Specifically, by exploring the high maneuverability of UAV, communication links can be established ubiquitously, especially in temporary hotspots, disaster areas, and complex terrains [13]. Due to the above benefits, numerous efforts have been endeavoured to the research of UAV assisted wireless communication design. For example, [14] proposed a suboptimal algorithm to solve the joint user association, UAV’s trajectory and power allocation design. [15] considered a solarpowered UAV system and a monotonic optimization based algorithm was developed to find the optimal UAV’s trajectory and power allocation. To benefit the advantages of multiantenna technique, [16] investigated the suboptimal beamforming design and UAV positioning to maximize the throughput of UAV assisted network. Considering the advantages of NOMA and UAV, the application of NOMA to the UAV assisted network was recently investigated in [17, 18]
. However, all the aforementioned works assumed that the channel state information (CSI) of the system can be perfectly acquired by the UAVs. In practice, due to, e.g., imperfect channel estimation and finite feedback, the UAVs can never have perfect CSI. Moreover, the wind would incur nonnegligible body jittering of the UAVs, which will also harm the acquisition of CSI
[19]. Therefore, it is of great importance to investigate the robust system design of UAV assisted wireless communication system under imperfect CSI assumptions. The robust beamforming design has been well studied in the orthogonal multiple access systems, e.g., [20, 21, 22], and recently, has been extended to the NOMA system, e.g., [23, 24, 25]. More recently, the robust beamforming design has also been introduced to the UAV enabled system, e.g., [26, 27]. Specifically, in [26], physical layer networking coding and multiple user detection approaches were proposed to combat the effects due to the absence of CSI, and [27] analyzed the power allocation problem in the singleantenna UAVNOMA systems.In this paper, we study the joint user clustering and robust beamforming design problem of a downlink multiuser multiUAV NOMA network with imperfect CSI assumption, where the users are spatially located by flowing the Poisson Cluster Process (PCP). Due to the integer variables relevant to user clustering, coupling effects of the downlink beamformers, and the infinitely many constraints bringing by the imperfect CSI, the resultant problem is challenging to solve. For computational complexity reduction, the original problem is divided into the user clustering subproblem and robust beamforming design subproblem. For the user clustering subproblem, unlike [10, 9, 28, 11]
, where traditional optimization approaches are used, we would like to resort to unsupervised learning based approaches to solve the user clustering problem by exploring the users’ position information. In particular, the users locates closely will be grouped into the same cluster. This is intuitively nontrival as the users locating together are likely to have similar channel conditions and thus can be efficiently served by a hovering UAV. As a first attempt, treating the users’ position information as feature data,
[29]proposed a kmeans based unsupervised clustering algorithm to solve the user clustering problem. However, the kmeans based approach is sensitive to the initial cluster centroids selection and improper initial cluster centroids may result in undesirable user clustering outcomes. Therefore, it is of great interests to design a more robust machine learning based approach to efficiently solve the user clustering problem. Under imperfect CSI assumption, the worsecase robust beamforming design problem, with signaltointerferencenoiseratio (SINR) and SIC constraints, is also investigated. Different to
[26, 27], which considered either the singlecell or singleantenna scenario, this work focus on a multiantenna scenario in the multicell interference channel. The coupling effects of the beamformers, both from inter and intracluster, make the formulated problem nonconvex and challenging to solve. Moreover, the imperfect CSI assumption make the considered problem much more complicated, due to the fact that each of the worstcase SINR or SIC constraints corresponds to an infinite number of nonconvex constraints. For more practical applications, we further investigate the algorithm to solve the robust beamforming design problem in a decentralized fashion. The contributions of this work are summarized as follows:
We formulate the joint user clustering and robust beamforming problem in a multiantenna UAVNOMA system under imperfect CSI assumptions. For computationefficiency consideration, the original problem is decoupled into two subproblems, i.e., the user clustering problem and the robust beamforming problem. To solve the user clustering problem efficiently, we propose to use a kmeans++ based unsupervised learning approach, which consists of a careful initial cluster centroids selection process and a standard k means based user clustering process.

To attain some useful insights on solving the robust beamforming design problem, we first consider a special case with perfect CSI. By equivalently transforming the associated problem into a secondorder cone programming, the relevant problem is solved optimally. Then, we focus on the problem in the general case with imperfect CSI. To simplify the corresponding problem, we first use the semidefinite relaxation (SDR) based method to transform the quadratic terms respect to beamformers into linear ones, and then, the Slemma is invoked to dealt with the infinitely many constraints caused by the imperfect CSI. By omitting the rankone constraints, the reformulated problem refers to the semidefinite programming (SDP) problem, which is convex and can be efficiently solved by the existing optimization tools. Further, to gain more insights on the proposed SDR based algorithm, a sufficient condition, under which the rankone optimality of the obtained solution can be guaranteed, is provided and the rankone optimality of the obtained solution is theoretically proved.

For practical applications, we investigate the decentralized approach to solve the robust beamforming design problem. By equivalent reformulations, the constraints set is decoupled into several independent subsets, each of which is solely related to a single UAV. Then, the alternating direction method of multipliers (ADMM) is applied to solve the reformulated problem efficiently.
Simulation results are presented to show the efficiency of the proposed transmission scheme and algorithms. The rest of the paper are organized as follow. In Section II, the considered system model and the transmission power minimization problem is introduced. Section III presents the kmeans++ based algorithm to solve the user clustering problem. In Section IV, the centralized robust beamforming design is investigated. While the decentralized robust beamforming design is studied in Section V. Simulation results and conclusions are given in Sections VI and VII respectively.
Notations
: Column vectors and matrices are denoted by boldfaced lowercase and uppercase letters, e.g.,
and ; , and stand for the sets of dimensional real and complex vectors and complex Hermitian matrices, respectively. denotes the identity matrix, and denotes an allzero (one) vector (matrix) with appropriate dimension. indicates that the element of are nonnegative. The superscripts , and describe the transpose, (Hermitian) conjugate transpose and pseudoinverse operations, respectively. and and represent the rank and trace of matrix, respectively. means that matrix is positive semidefinite (positive definite). denotes the Euclidean norm of vector . denotes the statistical expectation.Ii System Model and Problem Formulation
Consider a downlink communication network which consists of multiantenna UAVs and singleantenna users. Let and denote the index sets of UAVs and users, respectively. According to their locations, the users are grouped into nonoverlapping user clusters, i.e., where denotes the user set of the th cluster and is the number of users in the th cluster, satisfying . Each cluster is served by a hovering UAV. Without loss of generality, the users are assumed to locate in the same plane. Let and denote the coordinate vectors of user in the th cluster and UAV , respectively. In particular, is the flying height of UAV and , denote the coordinate of the centroid of the th cluster, which can be calculated by
(1) 
Iia User Location Model
Unlike the conventional assumption in the previous work that the users are uniformly distributed in the network
[30, 31], we focus on the scenario where the users’ locations are physically correlated in this work. This scenario characterizes the case where the system contains some hotspots, such as sports bar or lecture hall, etc [32]. In this case, the high maneuverability of UAV can be exploited to provide specific services for each hotspot. The Possion Cluster Process (PCP) is used to model the users’ distribution. Specifically, the PCP can be model by [33](2) 
where denotes the process of the parent points and denotes the offspring points process associated with the cluster center . We assume that the parent points are uniformly distributed in the network and the offspring points in the th cluster also follows the uniform distribution in a circular range, with radius of , around the cluster center.
IiB Channel Model
Note that the users in a cluster are severed by a UAV hovering above them, thus there is a high possibility that there exits a lineofsight (LoS) communication link between the UAV and its home user. Therefore, the channel between the UAV and its home user is modelled by the Rician fading channels. While, due to the blockages of the tallbuildings or trees, it is possible that no direct communication link exists between the UAV and its neighbouring users. Hence, the Rayleigh fading channel is more suitable for modelling the channels between UAV and its neighbouring users.
Let denote the channel from UAV to U. As discussed previously that the UAVs inevitably suffer from CSI errors in practice. Thus, imperfect CSI model is considered in this work. Let denote the preassumed CSI at the UAV for U. Then, the real CSI between UAV and U is given by
(3) 
where is the bounded CSI error associated with . In particular, the bounded CSI error can be modelled by
(4) 
where determines the range and shape of the CSI error. For instance, characterizes the popular spherical error model .
IiC Transmission Model and Problem Formulation
For spectral efficiency consideration, the NOMA protocol is applied to each user cluster. Follow the rational of NOMA, the signals of users in cluster are combined by using the superposition coding technique at UAV , then the users with stronger channel conditions will first remove the signals for the users with weaker channel conditions by invoking the SIC technique. Let denote the signal for user in cluster with . So, after superposition coding, the transmit signal of UAV is given by
(5) 
where is the beamformer for U. The received signal at U is give by
(6) 
where is the received additive white Gaussian noise at U
with zero mean and variance
. The first item in (6) denotes the desired signal of U, the second and third items in (6) denote the intracell and intercell interference, respectively.Without loss of generality, we assume that the users, in each cluster, are ordered by their channel gains in a descending manner, i.e., . Thus, according to the principle of NOMA, U would first remove the information, , for U for by using SIC, and then decoding its own information. Based on the above model, the SINR for decoding , , at U in the SIC process is given by
(7) 
After is removed from the received signal, the SINR at U for decoding is given by
(8) 
Based on the above model, the energyefficient joint userclustering and robust beamforming design problem can be formulated as
(9a)  
(9b)  
(9c)  
(9d)  
(9e) 
where is the transmission power budget of each UAV. (9b) represents that each user will be uniquely assigned into one cluster; (9c) guarantees the success of SIC procedure at each user and the QualityofService requirement for each user is given in (9d).
It is not difficult to verify that problem (9) is a nonconvex optimization problem due to the coupling of the quadratic beamforming vectors and also the channel uncertainty. More precisely, problem (9) is in fact an NPhard mixed integer programming problem and, thus, is unsolvable within polynomial time. To efficiently produce a high quality solution of problem (9), similar to [10, 9, 28, 11], we will decouple it into two subproblems, i.e., the user clustering problem and the robust beamforming design problem. However, unlike [10, 9, 28, 11], where the user paring problem is solved by the traditional optimization methods, in this work, we will use an unsupervised clustering based approach to efficiently produce a highquality clustering outcome with much lower computational complexity. In the next section, the details of the kmeans++ based unsupervised user clustering algorithm will be discussed.
Iii A KMeans++ based Approach for User Clustering
Notice that the users are spatially located in several clusters by following the PCP. It motivates us to design an efficient user clustering algorithm by exploring the distribution information of the users’ positions. As mentioned previously, to benefit the advantages of unsupervised learning and the users’ position information, [29] proposed a kmeans based algorithm to perform fast user clustering, but with the curse of the sensitivity to the initial centroids selection. In view of this, a new unsupervised learning based algorithm will be proposed in this work. To this end, by utilizing the users’ position information, we first rewrite the user clustering problem as the following Euclidean distanceoriented optimization problem
(10a)  
(10b)  
(10c)  
(10d) 
where collects the positions of all users, represents the centroids of the user clusters and indicates the cluster assignment of users. Specifically, implies that user is assigned to the th cluster, otherwise. Constraint (10b) indicates that each user can only be grouped into one cluster.
Notice that the constraint sets on and are decomposable. Thus, problem (10) can be efficiently solved by the alternating optimization (AO) based algorithm. In particular, by using the AO based algorithm, the user clustering outcome and the centroids of the clusters would be alternatively updated until the clustering converges. This is exactly the wellknown kmeans algorithm [34], which has been widely used to solve the data clustering problem, e.g., [29]. However, due to the nonconvex integer constraint in (10b), the kmeans algorithm is sensitive to the initial selection of the cluster centroids and may not always yield satisfactory clustering performance [35]. In view of this, the kmeans++ algorithm, which is an improvement of kmeans algorithm, will be used to solve the clustering problem in this work.
The ingredients of the kmeans++ method are twofolds: one is the initial cluster centroids selection process; the other is the standard kmeans method to find the final clustering outcome. The basic idea of the initial cluster centroids selection is as follow. The system first randomly selects a user as the first centroid. Secondly, the system needs to compute the distances from all the other users to this centroid, denoted by . Then, user will be selected as the second cluster centroid with probability . Thirdly, recompute the distances from all users to these two selected centroids. Then, let where and denote the distances from user to the first and second cluster centroid. Again, choose user as the third cluster centroid with probability . Repeat the above steps until all the centroids are selected. Finally, standard kmeans algorithm will be applied to solve the user clustering problem based on the chosen cluster centroids. The detailed steps of the kmeans++ based user clustering algorithm is summarized in Algorithm 1.
(11) 
From Algorithm 1, we can see that, one important trick of the initial cluster centroids selection is to choose the user that has a lager distance to the chosen centroids as the next centroid with a higher probability. This is intuitively reasonable, as larger distance between centroids results in a more robust user clustering outcome. Also, note that the kmeans++ algorithm chooses the next centroid with a probability, instead of choosing the user with the largest distance to the chosen centroids as the next centroid directly. This is to combat the effects from the noise point^{1}^{1}1The noise point in data clustering means that this point is far from all the other data points in the data set..
It is important to point out that, at the first glance, the computational complexity of the kmeans++ method would be higher than that of the standard kmeans method, as an additional initial cluster centroids selection process has been added. However, generally, the convergence behaviour of the kmeans++ can perform better than the kmeans method thanks to the careful selection of the initial cluster centroid selection. Although it is difficult to quantify theoretically, the performance efficacy, both in accuracy and speed, of the kmeans++ algorithm has been verified on various data sets, compared to the standard kmeans algorithm, see [35, section 6] for more details.
Iv Solve the Robust Beamforming Design Problem
Once the user clustering is determined, problem (9) boils down to a pure robust beamforming design problem. In this section, we will first optimally solve a special case of problem (9) under the perfect CSI assumption to gain some insights on the potential difficulty in solving problem (9). Then, we investigate the robust beamforming design present in problem (9) and an SDRbased algorithm is proposed to produce a suboptimal solution in an efficient way. Finally, the condition, under which the proposed SDRbased algorithm can generate an optimal solution, is studied.
Iva Optimal Design for the Case under Perfect CSI Assumption
Recall problem (9) under the assumption of perfect CSI as
(12a)  
(12b)  
(12c)  
(12d) 
By epigraph reformulation, problem (12) can be equivalently reformulated as
(13a)  
(13b)  
(13c)  
(13d)  
(13e) 
which is a secondorder cone programming and thus can be efficiently solved by the interior point based solver, e.g., CVX [36]. Now, one can realize that the challenge in solving the robust beamforming design problem is rather than the coupling effect of the quadratic beamformers, the channel uncertainty instead. In the next subsection, we will focus on the robust beamforming design problem.
IvB Suboptimal Design for the General Cases
In this section, we first propose to use the SDR method to relax the quadratic terms related to the beamforming vectors to linear ones. Then, we handle the obstacle bringing by the infinitely many SINR constraints. Finally, we also provide a sufficient condition under which the SDR will be tight.
B.1. SDRBased Suboptimal Design : To apply the SDR method, we first introduce a set of rankone matrix . Then, by ignoring the rankone constraint on the matrix, the robust beamforming problem can be relaxed as
(14a)  
(14b)  
(14c)  
(14d)  
(14e) 
which is a convex problem as the objective function and constraints are linear. However, it is still computationally intractable due to the infinite number of constraints. Next, to make problem (14) tractable, we propose the following lemma.
Lemma 1
Proof:
Based on Lemma 1, the SDP problem (14) can be equivalently reformulated as
(16a)  
(16b)  
(16c)  
(16d)  
(16e)  
(16f) 
which is an SDP and thus can be efficiently solved by CVX.
Remind that problem (14) is a relaxed version of the original robust beamforming design problem by ignoring the rankone constraints. Thus, one important issue in solving problem (14) is to verify whether the obtained matrices from solving problem (16
) are rankone. If it is true, then the optimal beamforming vectors can be obtained by simply applying singular value decomposition to the obtained matrices. Hence, it is interesting to explore the conditions under which solving problem (
16) can produce rankone solutions.B.2. RankOne Optimality Analysis: The following Lemma provides a condition that can guarantee the rank optimality of problem (16).
Lemma 2
Suppose that problem (16) is feasible, the rankone optimality can be guaranteed if for all , i.e., no intracell CSI error.
Proof:
The proof is relegated to Appendix B. ∎
V Decentralized Beamforming Design via ADMM
In the previous section, we propose to solve the robust beamforming design problem (16) by using SDR in a centralized fashion. However, solving problem (16), in this way, is based on the assumption that there is central controller to collect the entire CSI of the UAVs. In this case, the signalling overhead would become heavier as the increase of the number of users. Moreover, such a central controller may not be always available in practice. Therefore, it is of great importances to investigate the approaches that solve problem (16) in a decentralized way. In view of this, we propose to use the wellknown ADMM method to solve problem (16) in a decentralized fashion. However, applying ADMM method to problem (16) is not straightforward. Fortunately, by equivalent reformulations, problem (16) can be transformed to an appropriate structure under which the ADMM method can be applied. To this end, we first define the following auxiliary variables:
(17) 
where is the transmission power of UAV , and is the received intercell interference power from the other UAVs to U in the system. Then, problem (16) can be reformulated as
(18a)  
(18b)  
(18c)  
(18d)  
(18e)  
(18f)  
(18g) 
To apply ADMM, one important step is to decompose the constraint set into several independent subsets. Fortunately, we observe that, with problem unchanged, the subindices and in (18d) can be interchanged. Hence, the constraints from (18b) to (18f) can be decomposed into the following independent convex sets
(19) 
Further, define the following variables:
(20a)  
(20b) 
where and for . (20a) collects all the intercell interferences, while (20b) contains and (where ) that are relevant to UAV. Remind that , it can be verified that there exists a matrix , such that . Then, based on (V) and (20), problem (18) can be compactly rewritten as
(21a)  
(21b)  
(21c) 
Now, we are ready to apply ADMM. According to the rationale of ADMM, we solve the following augmented problem
(22a)  
(22b)  
(22c) 
where is a penalty parameter. As the first step of ADMM, we give the augmented Lagrangian of (22) as follows:
(23) 
where is the dual variable associated with constraint (22c). Then, according to the principle of ADMM, we have the following primal updates:
(24)  
(25) 
where is the iteration index. As it can be seen, problem (24) is convex, and thus, can be efficiently solved by CVX. While for quadratic optimization problem (25), we have the following closedform solution: