To support the exponentially growing data traffic and number of devices in future wireless network, non-orthogonal 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 , and the ultra-reliable and low-latency communications network . 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, 
studied the influences of user pairing on the outage probabilities of users in the NOMA system.
proposed to use a branch-and-bound based algorithm to solve the joint user pairing and power allocation optimally with the worst case computation complexity of NP-hard. To develop a computation-efficient algorithm for the user clustering problem, matching theory based heuristic algorithms were presented in[9, 10, 11]. Recently, to further reduce the computational complexity,  proposed to use the k-means 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 . 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 . Due to the above benefits, numerous efforts have been endeavoured to the research of UAV assisted wireless communication design. For example,  proposed a suboptimal algorithm to solve the joint user association, UAV’s trajectory and power allocation design.  considered a solar-powered 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 multi-antenna technique,  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 non-negligible body jittering of the UAVs, which will also harm the acquisition of CSI. 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 , physical layer networking coding and multiple user detection approaches were proposed to combat the effects due to the absence of CSI, and  analyzed the power allocation problem in the single-antenna UAV-NOMA systems.
In this paper, we study the joint user clustering and robust beamforming design problem of a downlink multi-user multi-UAV 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,
proposed a k-means based unsupervised clustering algorithm to solve the user clustering problem. However, the k-means 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 worse-case robust beamforming design problem, with signal-to-interference-noise-ratio (SINR) and SIC constraints, is also investigated. Different to[26, 27], which considered either the single-cell or single-antenna scenario, this work focus on a multi-antenna scenario in the multi-cell interference channel. The coupling effects of the beamformers, both from inter- and intra-cluster, 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 worst-case 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 multi-antenna UAV-NOMA system under imperfect CSI assumptions. For computation-efficiency 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 k-means++ 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 second-order 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 S-lemma is invoked to dealt with the infinitely many constraints caused by the imperfect CSI. By omitting the rank-one 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 rank-one optimality of the obtained solution can be guaranteed, is provided and the rank-one 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 k-means++ 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.
: 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 all-zero (one) vector (matrix) with appropriate dimension. indicates that the element of are non-negative. The superscripts , and describe the transpose, (Hermitian) conjugate transpose and pseudo-inverse 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 multi-antenna UAVs and single-antenna users. Let and denote the index sets of UAVs and users, respectively. According to their locations, the users are grouped into non-overlapping 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
Ii-a 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 . 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 
where denotes the process of the parent points and denotes the off-spring points process associated with the cluster center . We assume that the parent points are uniformly distributed in the network and the off-spring points in the -th cluster also follows the uniform distribution in a circular range, with radius of , around the cluster center.
Ii-B 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 line-of-sight (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 tall-buildings 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 pre-assumed CSI at the UAV for U. Then, the real CSI between UAV and U is given by
where is the bounded CSI error associated with . In particular, the bounded CSI error can be modelled by
where determines the range and shape of the CSI error. For instance, characterizes the popular spherical error model .
Ii-C 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
where is the beamformer for U. The received signal at U is give by
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 intra-cell and inter-cell 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
After is removed from the received signal, the SINR at U for decoding is given by
Based on the above model, the energy-efficient joint user-clustering and robust beamforming design problem can be formulated as
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 Quality-of-Service 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 NP-hard 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 high-quality clustering outcome with much lower computational complexity. In the next section, the details of the k-means++ based unsupervised user clustering algorithm will be discussed.
Iii A K-Means++ 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,  proposed a k-means 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 distance-oriented optimization problem
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 well-known k-means algorithm , which has been widely used to solve the data clustering problem, e.g., . However, due to the non-convex integer constraint in (10b), the k-means algorithm is sensitive to the initial selection of the cluster centroids and may not always yield satisfactory clustering performance . In view of this, the k-means++ algorithm, which is an improvement of k-means algorithm, will be used to solve the clustering problem in this work.
The ingredients of the k-means++ method are two-folds: one is the initial cluster centroids selection process; the other is the standard k-means 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 k-means algorithm will be applied to solve the user clustering problem based on the chosen cluster centroids. The detailed steps of the k-means++ based user clustering algorithm is summarized in Algorithm 1.
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 k-means++ 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 point111The 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 k-means++ method would be higher than that of the standard k-means method, as an additional initial cluster centroids selection process has been added. However, generally, the convergence behaviour of the k-means++ can perform better than the k-means 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 k-means++ algorithm has been verified on various data sets, compared to the standard k-means 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 SDR-based algorithm is proposed to produce a suboptimal solution in an efficient way. Finally, the condition, under which the proposed SDR-based algorithm can generate an optimal solution, is studied.
Iv-a Optimal Design for the Case under Perfect CSI Assumption
Recall problem (9) under the assumption of perfect CSI as
By epigraph reformulation, problem (12) can be equivalently reformulated as
which is a second-order cone programming and thus can be efficiently solved by the interior point based solver, e.g., CVX . 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.
Iv-B 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. SDR-Based Suboptimal Design : To apply the SDR method, we first introduce a set of rank-one matrix . Then, by ignoring the rank-one constraint on the matrix, the robust beamforming problem can be relaxed as
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.
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 rank-one constraints. Thus, one important issue in solving problem (14) is to verify whether the obtained matrices from solving problem (16
) are rank-one. 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 rank-one solutions.
B.2. Rank-One Optimality Analysis: The following Lemma provides a condition that can guarantee the rank optimality of problem (16).
Suppose that problem (16) is feasible, the rank-one optimality can be guaranteed if for all , i.e., no intra-cell CSI error.
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 well-known 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:
where is the transmission power of UAV , and is the received inter-cell interference power from the other UAVs to U in the system. Then, problem (16) can be reformulated as
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
Further, define the following variables:
where and for . (20a) collects all the inter-cell 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
Now, we are ready to apply ADMM. According to the rationale of ADMM, we solve the following augmented problem
where is a penalty parameter. As the first step of ADMM, we give the augmented Lagrangian of (22) as follows:
where is the dual variable associated with constraint (22c). Then, according to the principle of ADMM, we have the following primal updates: