Geometrical Synthesis of Sparse Antenna Arrays Using Compressive Sensing for 5G IoT Applications

One of the main targets of the forthcoming fifth-generation (5G) cellular network will be the support of the communications for billions of sensors and actuators, so as to finally realize the Internet of things (IoT) paradigm. This pervasive scenario unavoidably requires the design of cheap antenna systems with beamforming capabilities for compensating the strong attenuations that characterize the millimeter-wave (mmWave) channel. To address this issue, this paper proposes an iterative algorithm for sparse antenna arrays that enables to derive the number of elements, their amplitudes, phases, and positions in the presence of constraints on the far-field pattern. The algorithm, which relies on the compressive sensing approach, is formulated by transforming the original nonconvex optimization problem into a convex one. To prove the suitability of the conceived solution for 5G IoT mmWave applications, numerical examples and comparisons with other existing methods are provided, considering synthesis problems with different pattern and aperture specifications.


Introduction
The idea of combining the radiation of multiple antennas to obtain specific pattern shape dates back to the beginning of the 19th century, when the history of wireless communications began [1,2]. Nowadays, antenna arrays are widely used in a large variety of scenarios involving far-and near-field focusing/multifocusing applications for satellite, cellular, vehicular, and sensor networks [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17]. However, even if the array technology is already widespread, its importance is expected to further increase, since the directionality of the communications will represent a basic enabling functionality of the forthcoming fifth-generation (5G) and Internet of things (IoT) systems [18][19][20][21]. This forecast is motivated, on one hand, by the expected presence of a huge number of active devices (smartphones, sensors, actuators), and, on the other hand, by the adoption of the millimeter-wave (mmWave) bands. In fact, the need of attenuating the reciprocal interference among the 5G equipments and of compensating the significant attenuations that characterize the mmWave channel implies the implementation of multi-antenna systems satisfying compactness and performance constraints. These constraints must be of course combined with the possible reduction in the unit price of a device, so as to better match the market demand. An initial way to achieve this requirement is to minimize the number of radiators in each antenna array. However, due to the so crowded electromagnetic environment, this task cannot be accomplished at the expense of the desired radiation characteristics. Therefore, the array synthesis problem addressed by a 5G antenna designer is very hard to solve, since it requires to optimize not only the excitations of the array elements, but also their positions and, possibly, their number. This situation identifies the typical problem addressed in the research field represented by the design of sparse antenna arrays [22]. This latter issue has attracted the interests of the research community since the 1960s [23][24][25][26][27], regaining attention in the more recent years thanks to its applicability to the emerging communication scenarios. Accordingly, several methods have been proposed in the literature for the synthesis of sparse antenna arrays by considering both stochastic and deterministic approaches [28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43]. Stochastic methods can rely on different strategies, including genetic algorithms [28,29], particle swarm optimization [30], differential evolution [31], and nature inspired techniques, such as the ant [32], whale [33], and grey wolf [34] optimization algorithms. Also deterministic methods have been developed by moving from different strategies. In particular, a noniterative algorithm is proposed in [35], where the matrix pencil method is adopted to approximate a desired pattern using a nonuniform linear antenna array with a reduced number of elements. In [36], the same authors present a deterministic three-step procedure able to optimize the magnitudes, phases, and locations of the elements of a linear antenna array radiating a shaped power pattern. In [37], the alternating projection algorithm is applied in conjunction with a special technique, which iteratively places the radiating elements with the aim of synthesizing a sparse array starting from a power pattern mask specification. In [38], a low-complexity deterministic method for the synthesis of phase-only beam scanning linear aperiodic arrays is proposed, which optimizes the positions, the amplitudes, and the phases of the array elements. Convex optimization techniques are instead applied to the synthesis of sparse antenna arrays in [39][40][41][42][43], where, the compressive sensing (CS) strategy is employed. The CS technique was originally introduced in the field of image processing and signal reconstruction [44], but nowadays its applications have hugely spanned. Thus, many different types of engineering problems have been solved by the CS approach, some of which specifically belongs to the electromagnetics' research context, such as the diagnosis and synthesis of antenna arrays, the estimation of the direction of arrivals, and the solution of inverse scattering and radar imaging problems [45]. With particular reference to the synthesis of sparse arrays, the main advantage of CS consists in the possibility to convert a nonconvex formulation into a convex one, thus facing a complex problem in an acceptable processing time.
According to the above observations, this paper considers the synthesis of sparse antenna arrays for mmWave applications by developing an iterative algorithm based on the CS optimization. The algorithm is conceived to address the general case in which a fixed grid structure is not available, and so nonlinear problems characterized by complex formulations have to be managed, even when simple array geometries and not challenging radiation constraints are involved. In fact, the aim is to enable the 5G antenna designer to synthesize not only the array excitations (in amplitude and phase), but also the number and the location of the elements in the presence of requirements on the far-field pattern.
In this sense, the main advantage of the presented method with respect to the previously proposed CS-based ones, consists in the possibility to impose that the synthesized pattern belongs to a predefined mask. Differently from the conventional pattern matching approach, where a specific shape is imposed both inside the main lobe (region of interest) and outside it (region not of interest), the here proposed algorithm applies in detail the shape constraint to the sole main lobe, imposing just an upper limit to the other angular regions. Moreover, the main lobe shaping is realized independently of the pattern phase. Thus, a power pattern problem is here addressed instead of a field synthesis one analyzed in [40]. In this way, the available degrees of freedom are better exploited. Besides, the pattern requirements are matched by reformulating the initial nonconvex problem into a convex one in order to enable its mathematical tractability in an acceptable CPU time. To prove the effectiveness of the presented method, numerical examples adherent to the 5G IoT context and comparisons with existing approaches are presented and discussed.
The remaining of the paper is organized as follows. Section 2 formulates the addressed problem. Section 3 presents the developed algorithm. Section 4 discusses the numerical results. Section 5 summarizes the most relevant conclusions.
Notation. Throughout the paper the following notation is used: (·) T denotes the transpose operator, (·) * denotes the complex conjugate, j denotes the imaginary unit, and ∠x denotes the argument of x.

Problem Formulation
With reference to a Cartesian system O(x, y, z), consider an antenna array of arbitrary geometry composed by a generic (and possibly large) number of elements N, which are located at the positions r n = x nx + y nŷ + z nẑ , for n = 1, . . . , N, wherex,ŷ andẑ represent the unit vectors of the corresponding Cartesian axes. As usual, in spherical coordinates, denote by θ the angle measured from the z-axis and by φ the azimuth angle. The radiation pattern produced by this arbitrary array at the generic direction r = sin θ cos φx + sin θ sin φŷ + cos θẑ is given by: where w = [w 1 , . . . , w N ] T represents the column vector of the complex excitations of the array elements, f n (r) is the element pattern of the n-th array element, and k = 2π/λ is the wavenumber, being λ the carrier wavelength.
Concerning (1), two main aspects should be taken into account. Firstly, the geometry of the array, which may be linear, planar or even conformal, is usually specified by some shape/size constraints in the physical design. When the CS strategy is adopted, the application of these constraints allows one to identify the possible positions of the N elements as the candidate positions, while the final array will be composed by a very reduced number of elements, suitably chosen among these N candidates. The second aspect that should be considered is that also the pattern generated by (1) is always required to satisfy specific shape constraints, which can be properly modeled by a suitably defined mask. Accordingly, the constraints on the pattern can, in principle, be expressed as: where M low (r) and M up (r) are real positive functions representing, respectively, the lower and upper bound of the mask. The problem addressed in this paper is that of selecting, among the N candidate positions, the lowest number of array elements, their positions and excitations in such a way as to obtain a radiation pattern compliant with (2). By consequence, the here addressed problem can be mathematically formulated in compact form as follows: where w 0 denotes the l 0 -norm of w. More precisely, w 0 counts the nonzero components of w, corresponding to the elements that will result physically active at the end of the synthesis process.

Synthesis Method
As it can be immediately observed, the problem in (3) is in general nonlinear and nonconvex, thus extremely hard to solve. The CS approach is specifically useful for this kind of situations, but requires a proper reformulation of (3). To this aim, consider an iterative procedure, which, at the generic k-th iteration, solves a minimization problem whose objective function is given by the weighted l 1 -norm (or weighted Manhattan distance) [39]: in which, for n = 1, . . . , N, k = 1, 2, . . . , the weights: are introduced to improve the final result by properly selecting a suitable parameter ε. The objective function in (4) along with the weight definition in (5) replaces the original objective function in (3a), and is known to produce sparse solutions [46]. However, this novel problem remains difficult to solve because of the nonlinearity and nonconvexity of (2)-(3b). To address this second aspect, let outline some observations regarding the radiation pattern requirements that are commonly imposed. Usually, a desired pattern is characterized by: (i) a main beam region (MBR) with a shape that must be exactly matched, (ii) a maximum allowed level for the sidelobe region (SLR), and (iii) a (possible) null region (NR). For these two latter requirements, the lower bound of the mask is not of concern, since the lower is the pattern, the more appreciable is the result. Therefore, according to [40], a suitable strategy to handle the constraint in (2)-(3b) is to separately impose the requirement on the MBR and those on the SLR and NR. Accordingly, (2) is subdivided into the following two constraints: where F d (r) represents the desired main lobe function [40]. Now, it is worth to note that the usage of (6a) leads to a field synthesis problem and not to a power synthesis one. This implies that the available degrees of freedom are in part wasted to approximate the array phase pattern, which is usually not of interest. To overcome this issue, one can first recall the following property holding for two arbitrary complex numbers z 1 and z 2 : This property suggests that (6a) is equivalent to a power synthesis problem if and only if the phase of the array pattern equals the phase of the desired pattern. Hence, by considering (4)- (7), the original problem in (3) can be iteratively solved by minimizing, at the k-th iteration, the constrained function: where δ(r) is a real function defining the accuracy required in the MBR at each direction, and: is the desired MBR function updated, at the k-th iteration, according to the phase of the pattern synthesized at the (k − 1)-th iteration and to a suitable real function F 0 d (r) defining the MBR shape. Notably, this latter formulation represents a second-order cone program (SOCP) problem, which can be solved by the use of freely available software routines, as, for example, the Matlab-based CVX [47]. It is also worth to observe that, although, at each iteration, (8a)-(8c) formally belongs to the class of field synthesis problems, the iterative modification in (9) of the constraint in (8b) actually leads to a power synthesis problem. Importantly, note that it is this iterative modification of the MBR constraint that allows one to better exploit the degrees of freedom of the problem and thus to improve the final results. Moreover, to refer all requirements in terms of mask specifications, one can define the functions identifying the accuracy and the MBR shape, respectively, as: This operation enables to finally identify the iterative algorithm that solves the original synthesis problem formulated in (3). The development of the algorithm is detailed in Table 1.
Step 1 specifies the problem's constraints (candidate positions, bounds of the mask, accuracy, SLR, NR, and initial MBR shape).
Step 2 initializes the iteration and the excitations.
Step 3 updates the iteration, the pattern, the weights, and the phase of the MBR shape. This latter update constitutes the major novelty of the proposed approach.
Step 4 solves the SOCP problem at the present iteration.
Step 5 identifies the stop condition, which requires that the number of nonzero elements of the excitation vector does not change in three consecutive iterations. Alternatively, the stop condition might also be formulated by selecting a maximum pre-defined number of iterations. Of course, during the evolution of the algorithm, none element of w k does exactly vanish, but the generic n-th element of the array is considered as zero, that is, is assumed absent, if the amplitude of its excitation |w k n | is lower than the threshold ε. Step 1 Problem specifications (i) Initial set of candidate positionsr n (n = 1, . . . , N); (ii) Upper M up (r) and lower M low (r) bounds of the mask; (iii) MBR, SLR, and NR; (iv) Accuracy δ(r) and MBR shape F 0 d (r) using (10a) and (10b), respectively.

Stop condition
If k ≥ 3 and w k 0 = w k−1 0 = w k−2 0 w k identifies the elements of the final sparse array and their excitations; Exit else Go to Step 3 end

Results and Discussion
Five numerical examples are provided to prove the effectiveness of the proposed algorithm, by considering both two-dimensional (2D) problems (i.e., azimuth synthesis only) and three-dimensional (3D) ones (i.e., zenith-azimuth synthesis). Four of these examples are taken from the literature, in order to have a direct benchmark with the state-of-the-art methods in terms of element saving, while the fifth example is specifically conceived to include a more evolved multi-ring array geometry. Note that the first four problems are intentionally selected from examples already proposed in previous papers, since, according to the approach commonly adopted in the array synthesis research field, a direct comparison between the proposed and the existing algorithms can be immediately carried out. Furthermore, to put into evidence the relevance of all the solved problems for the IoT scenario, each example is described including specific references in which the employed array is exploited for 5G sensor applications. The algorithm is implemented using Matlab R2018b and CVX on a personal laptop having 8 GB RAM, and all results are obtained adopting a threshold ε = 10 −2 [40].  Table 2.

First Example
The first numerical example considers the 2D problem addressed in [39], which involves a linear array of isotropic elements lying on the z-axis and characterized by an aperture equal to 20λ. This array is required to radiate in the zenith domain a flat-top pattern, in which the MBR is given by the set Θ MBR = {θ : 70 • ≤ θ ≤ 110 • } with an allowed ripple ρ t equal to 0.4455 dB, while the SLR is defined by the set Θ SLR = {θ : 0 • ≤ θ ≤ 65 • ∪ 115 • ≤ θ ≤ 180 • } with a maximum sidelobe level equal to −30 dB. The initial set of candidate positions is selected by regularly spacing the isotropic radiators on the available aperture with an inter-element distance of λ/100, resulting in a starting value N = 2001. The linear array geometry has been widely applied to IoT scenarios [19,20,48], thanks to its simplicity of realization and conformability to objects having the length as the prevailing dimension.
The pattern obtained using the proposed algorithm is reported in Figure 1. Table 2, instead, lists the element positions normalized with respect to λ and the excitations normalized with respect to w 10 . From these values one may first observe that, being the desired pattern symmetrical, the active elements are symmetrically displaced and the corresponding excitations are real numbers. Of course, the central element n = 10 lies in position z 10 = 0 and has a normalized excitation equal to one. Interestingly, among the N = 2001 candidates, just 19 elements have been sufficient to match the pattern requirements. This represents a considerable reduction of the final number of elements as compared to [39], in which the final array was composed by 31 elements, a result already better than the 41 elements calculated in [49,50], where this problem was originally developed.   Table 3.
A possible further reduction of the number of elements might be achieved by imposing a further minimum inter-element distance control, which in this case can be considered suitable. In fact, one may notice from Table 2 that the elements n = 8 and n = 9, and, similarly, the elements n = 11 and n = 12, are very close to each other. In particular, their distance is exactly equal to the grid step λ/100, thus making difficult the practical realization of the array. The proposed refinement consists in replacing the pairs (8,9) and (11,12) with two single elements lying in the middle points (z 8 +z 9 )/2 = −0.645λ and (z 11 +z 12 )/2 = 0.645λ, and selecting, for the respective excitations, the values w 8 +w 9 = 0.6351 and w 11 +w 12 = 0.6341. In this way, one can obtain an array of only 17 elements that radiates the pattern in Figure 1 identified by the red dashed line, which is characterized by a very limited degradation with respect to the originally synthesized one (blue solid line). As a final observation, it is worth to note that the here derived original pattern has been obtained after six iterations that have required 254 s. Thus, just a few minutes have been sufficient to achieve the presented result.

Second Example
The second example is still a 2D flat-top synthesis problem taken from [39]. The problem considers the same linear array and candidate elements of the first example, but, differently, the array is composed by directive radiators, which are more likely to be mounted on actual IoT sensors [19,20,48]. More precisely, the single-element pattern is modeled by the function f n (θ) = sin θ, representing a short dipole parallel to the z-axis. Moreover, the desired MBR is no more symmetrical with respect to the broadside direction, but is defined by the shifted set Θ MBR = {θ : 50 • ≤ θ ≤ 90 • } with an allowed ripple ρ t equal to 1 dB, while the SLR is defined by the set The pattern derived using the developed method is shown in Figure 2, while Table 3 reports the normalized element positions and the complex excitations. In this second example, the desired pattern is not symmetrical, and hence the positions are also not symmetrical and the excitations are not real. The provided results confirm the satisfactory behavior of the proposed technique, since the final array, derived after 13 iterations that have required a CPU time of 548 s (less than 10 min), consists of just 18 active elements. A significant improvement as compared to the 25 elements obtained in [39]. Beside the element reduction, this example proves that the presented method is suitable not only when broadside patterns and isotropic sources are assumed, but also when more general scenarios, characterized by steered main beams and directive radiators, have to be realistically managed.   Table 4.

Third Example
The third example involves the 2D synthesis problem originally developed in [51] that adopts a linear array of aperture 7.5λ composed by isotropic elements. In this case, the algorithm is required to generate a not symmetrical cosecant-like pattern in the MBR defined by the set Θ MBR = {θ : 98 • ≤ θ ≤ 135 • }, and to impose a NR given by the set Θ NR = {θ : 66 • ≤ θ ≤ 88 • } in which no more than −30 dB of radiation are allowed. Besides, the SLR is defined by the set Θ SLR = {θ : 0 • ≤ θ ≤ 66 • ∪ 143 • ≤ θ ≤ 180 • } with a maximum sidelobe level equal to −20 dB, and the pattern obtained in [36] for the same example is used as the desired pattern in (9) in the MBR. Also in this case, the candidate positions are chosen by regularly spacing the radiators on the available aperture with an inter-element distance equal to λ/100, resulting now in a starting value N = 751. Note that, with reference to forthcoming pervasive communication networks, the cosecant-like pattern can be of specific interest for IoT sensors, since it enables to radiate a constant power on a given angular region, thus realizing a uniform covering of that region for monitoring applications [52,53].
The pattern calculated employing the presented approach is reported in Figure 3, while Table 4 shows the normalized element positions and the complex excitations. These results reveal the capability of the developed algorithm to strictly approximate the desired pattern, even in the presence of a wide null region and of a main beam requiring a specific not-flat shape. However, the main advantage of the designed method remains the low number of resulting active elements as compared to the previously proposed techniques. In fact, the original solution in [51] was characterized by 16 elements, and those in [36] and [40] by 13 elements, while the here conceived algorithm has provided, after just four iterations, a final array consisting of only 12 active elements. Furthermore, the time necessary to achieve this result has been approximately equal to 55 s, thus lower than a minute.   Table 5.

Fourth Example
The fourth example is a 3D zenith-azimuth synthesis problem originally developed in [39] that involves a square array of isotropic radiators lying in the x − y plane and having side equal to 5λ. The mask specifications are given in terms of the variables u = sin θ cos φ and v = sin θ sin φ. In particular, the MBR, in which the allowed ripple ρ t is equal 1 dB, is defined by the 2D set Ω MBR = {(u, v) : u 2 + v 2 ≤ 1/25} while the SLR, in which the maximum allowed level is equal to −25.85 dB, is identified by the 2D set Ω SLR = {(u, v) : u 2 + v 2 ≥ 4/25}. The initial set of the candidate positions is composed by a regular grid of elements spaced by λ/4 in both directions and covering the available aperture. This leads to a starting value N = 441. This fourth example is of interest both for 5G base stations (BSs) and gigabit-WiFi access points (APs) [54], which are, on one hand, characterized by a planar structure that enables the possibility to host many radiating elements [55][56][57], and, on the other hand, have to manage the problem of initial user access, which requires a wide main beam to avoid too long searching procedures for identifying the region of space where a given user lies [58].  The contour plot of the pattern derived using the proposed algorithm is shown in Figure 4, while Figure 5 and Table 5 report the initial grid (red cross) with the finally active elements (blue circles) and the real excitations, respectively. For the correspondence between Figure 5 and Table 5, the active elements have been numbered from the bottom to the top and from the left to the right. For readability reasons, Figure 5 illustrates the sole four elements that enable the reader to infer the order. From these results, one may notice that 60 elements have been obtained. A considerable reduction as compared to the 76 derived in [39]. The computational time required to achieve the here provided solution is equal to 2032 s, corresponding to 13 iterations.  Figure 7 and the excitations are listed in Table 6.

Fifth Example
The fifth and last example is a 3D zenith-azimuth synthesis problem developed adopting a planar circular array with multiple rings radiating a very narrow beam. The multi-ring array, having a maximum radius equal to 10λ, consists of isotropic radiators lying in the x − y plane. The initial set of the candidate positions is composed by N = 1308 elements regularly spaced on 20 rings having center in the origin of the reference system and radii R i = 0.5(1 + i)λ, i = 0, . . . , 19, in such a way as d min ≥ λ/2. The mask specifications, given in terms of the variables u and v, have the MBR defined by the 2D set Ω MBR = {(u, v) : u 2 + v 2 ≤ 1/400} and the SLR identified by the 2D set Ω SLR = {(u, v) : u 2 + v 2 ≥ 1/100}. The allowed ripple in the MBR is ρ t = 1 dB and the maximum allowed level in the SLR is equal to −15 dB. This example refers to a scenario in which the 5G BS or the gigabit-WiFi AP (both suitable to host planar array structures), has already accomplished the initial access with the user and intends to generate a highly directional communication for realizing a high-capacity link [38,59,60].
The contour plot of the pattern derived using the proposed algorithm is shown in Figure 6, while Figure 7 and Table 6 report the initial grid (red cross) with the finally active elements (blue circles) and the real excitations, respectively. For the correspondence between Figure 7 and Table 6, the active elements have been numbered starting from the x−axis counterclockwise from the outermost to the innermost ring. For readability reasons, Figure 7 illustrates the numbers of the four elements that enable the reader to infer the order. The computational time required to achieve the here provided solution is higher with respect to the previous example (14,691 s, corresponding to eight iterations). This is due to the finer grid required to sample a so narrow beam in the u − v domain. This last example and the previous one prove that also 3D synthesis problems may be successfully dealt with by the conceived CS-based approach, considering both flat-top and pencil beam pattern requirements.
A final comment regarding the presented results concerns the mutual coupling effects, which, in the research field covered by the geometrical synthesis of antenna arrays, must be considered once the final positions have been estimated. In fact, the initial set of positions does not represent real elements, but only candidate ones. Hence, when the selected elements change, also the coupling effects change. Consequently, one should evaluate the coupling at each iteration. Clearly, this is not a practical approach and hence the coupling effects may be considered at the end of the procedure. In this case, if the pattern distortion is not acceptable, the element excitations may be modified with any suitable algorithm for the synthesis of conformal arrays, as for example [49,50,[61][62][63][64].

Conclusions
A CS-based iterative procedure for the power synthesis of sparse arrays has been presented. The proposed algorithm relies on the solution of a sequence of SOCP problems with the aim of minimizing the number of radiators of an array composed by a large number of candidate elements. The constraints of the minimization problem have been formulated in such a way as to be convex, and to approximate a power pattern synthesis problem in the entire visible region. The constraints formulation, and, in particular, their iterative modification constitutes the original contribution, which allows one to improve the performance of the previously proposed CS-strategies for sparse array applications. Different numerical examples involving linear and planar structures have been discussed, obtaining, in all cases, the matching of the pattern requirements, and, for the cases involving the comparison with the previous existing solutions, a lower final number of active elements. Funding: This research was partly funded by the Italian Ministry of University and Research (MIUR) within the project FRA 2018 (University of Trieste, Italy), entitled "UBER-5G: Cubesat 5G networks -Access layer analysis and antenna system development."

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript: