*Article* **Joint User Scheduling and Hybrid Beamforming Design for Massive MIMO LEO Satellite Multigroup Multicast Communication Systems**

**Yang Liu <sup>1</sup> , Changqing Li 2,\*, Jiong Li <sup>2</sup> and Lu Feng <sup>1</sup>**


**Abstract:** In the satellite multigroup multicast communication systems based on the DVB-S2X standard, due to the limitation of the DVB-S2X frame structure, user scheduling and beamforming design have become the focus of academic research. In this work, we take the massive multi-input multi-output (MIMO) low earth orbit (LEO) satellite communication system adopting the DVB-S2X standard as the research scenario, and the LEO satellite adopts a uniform planar array (UPA) based on the fully connected hybrid structure. We focus on the coupling design of user scheduling and beamforming; meanwhile, the scheme design takes the influence of residual Doppler shift and phase disturbance on channel errors into account. Under the constraints of total transmission power and quality of service (QoS), we study the robust joint user scheduling and hybrid beamforming design aimed at maximizing the energy efficiency (EE). For this problem, we first adopt the hierarchical clustering algorithm to group users. Then, the semidefinite programming (SDP) algorithm and the concave convex process (CCCP) framework are applied to tackle the optimization of user scheduling and hybrid beamforming design. To handle the rank-one matrix constraint, the penalty iteration algorithm is proposed. To balance the performance and complexity of the algorithm, the user preselected step is added before joint design. Finally, to obtain the digital beamforming matrix and the analog beamforming matrix in a hybrid beamformer, the alternative optimization algorithm based on the majorization-minimization framework (MM-AltOpt) is proposed. Numerical simulation results show that the EE of the proposed joint user scheduling and beamforming design algorithm is higher than that of the traditional decoupling design algorithms.

**Keywords:** LEO satellite communications; massive MIMO; multigroup multicast; user scheduling; hybrid beamforming; robust; joint design; energy efficiency

## **1. Introduction**

In recent years, LEO satellite communication systems have played an increasingly important role in wireless communication networks [1]. However, facing the high-performance requirements of future wireless communication systems, such as higher spectrum efficiency (SE) and increased EE, the performance of LEO satellite communication systems needs to be improved [2]. Applying the massive MIMO technology to LEO satellites is a good choice, and with the advantage of 5G technology [3] and the spatial multiplexing principle of MIMO technology [4], the performance of LEO satellite communication systems can be further improved. Meanwhile, by using the high-precision multiple beams generated by the massive MIMO technology and aggressive full frequency reuse scheme among beams, the performance of the communication system can be greatly improved. However, the full frequency reuse scheme will cause severe inter-beam interference [5], and adopting the beamforming design at the LEO transmitter side can efficiently manage it. In addition, the super-frame structure of multibeam satellite communications standards such as DVB-S2X [6] needs to apply the same beamformer to multiple users that share the same frame.

**Citation:** Liu, Y.; Li, C.; Li, J.; Feng, L. Joint User Scheduling and Hybrid Beamforming Design for Massive MIMO LEO Satellite Multigroup Multicast Communication Systems. *Sensors* **2022**, *22*, 6858. https:// doi.org/10.3390/s22186858

Academic Editor: Josip Lorincz

Received: 10 August 2022 Accepted: 7 September 2022 Published: 10 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Therefore, the multigroup multicasting principle can be used in the beamformer design. In this paper, we concentrate on the massive MIMO LEO satellite communication system forward link multigroup multicasting beamforming scheme [7].

In the beamforming design of the massive MIMO LEO satellite multigroup multicast communication system, the following issues need to be considered:


#### **2. Related Works and Main Contributions**

#### *2.1. Related Works*

There is extensive literature regarding user scheduling and beamforming design in the wireless multigroup multicast communication system. In Ref. [11], based on the perfect CSI, taking the throughput maximization as the optimization objective the authors studied the precoding design of the multibeam satellite communication system, and proposed a decoupling scheme of the user scheduling and beamforming design. In Ref. [12], considering the influence of CSI errors and taking the minimizing transmission power as the optimization objective, the robust multigroup multicast transmission scheme of the multibeam satellite communication system was investigated, and the low complexity beamforming algorithm and the user grouping algorithm were proposed, but the length limit of the DVB-S2X frame was not considered. In Ref. [13], based on the perfect CSI and taking the maximizing SE as the optimization objective, the multigroup multicast transmission design scheme of the frame-based multibeam satellite communication system was investigated, and the authors proposed a joint design scheme of the user scheduling and beamforming. However, the influence of CSI errors was not considered, and the user grouping algorithm was simple. In Ref. [14], based on the perfect CSI, the authors studied the user scheduling problem of the multicast transmission in the high-throughput satellite communication system. The user scheduling was decoupled into intra-beam scheduling and inter-beam scheduling, and the correlation degree was calculated by using the equivalent CSI; therefore, the interaction of

intra-beam and inter-beam scheduling cannot be fully considered. In Ref. [15], the authors studied the user scheduling problem of the multibeam satellite communication system but did not consider the inter-beam interference caused by scheduling. In Ref. [16], based on the perfect CSI, the authors studied the joint scheduling and beamforming design problem for multiuser MISO downlink, and the message-based user grouping and scheduling algorithm was mainly proposed, but the impact of user grouping on system performance was not fully considered. In Ref. [17], the user scheduling and hybrid beamforming design of the massive MIMO orthogonal frequency division multiple access (OFDMA) communication system was studied, in scheme design. First, a joint design algorithm of user scheduling and analog beamforming was proposed; then, the digital beamforming matrix was solved by the weighted minimum mean-square error (WMMSE) algorithm. In Ref. [18], the authors studied the design of user scheduling and subcarrier allocation in the downlink of a massive MIMO OFDMA communication system and proposed a hybrid beamforming scheme. First, based on the optimal solution of digital beamforming, the analog beamforming matrix was obtained by a singular value decomposition algorithm. Then, the authors proposed an algorithm to solve the digital beamforming matrix and its corresponding scheduling users.

Most of the above research took the geosynchronous earth orbit (GEO) satellite or the terrestrial cellular network as the research object, and less used the LEO satellite. In addition, the optimization objectives mainly focused on maximizing SE and minimizing transmission power, and less on the EE optimization. Meanwhile, the analysis of CSI errors was insufficient, as some research considered the CSI errors, but did not analyze the influence of the Doppler shift. In terms of the user scheduling, the common idea was to adopt the decoupling scheme of user scheduling and beamforming design, without fully considering the coupling relationship between the user scheduling and the beamforming design.

#### *2.2. Main Contributions*

Inspired by the above research, we focus on the downlink transmission design of the massive MIMO LEO satellite multigroup multicast communication system. In scheme design, comprehensively considering the influence of CSI errors caused by the residual Doppler shift and the phase disturbance, we mainly investigate the robust joint user scheduling and hybrid beamforming design to maximize the system EE. Meanwhile, we take the constraints of the transmission power and QoS into account. The main works are summarized as follows:


• For the solution of the digital beamforming matrix and the analog beamforming matrix in the hybrid beamformer, the MM-AltOpt algorithm is proposed. trix in the hybrid beamformer, the MM-AltOpt algorithm is proposed. **3. System Model and Problem Formulation**

• For the solution of the digital beamforming matrix and the analog beamforming ma-

• For the DC programming problem, an iterative optimization algorithm based on the CCCP framework is proposed. For the rank-one matrix constraint introduced by the

formed into a difference of convex (DC) programming problem.

SDP algorithm, a penalty iterative algorithm is adopted.

form into a concave function, and some nonconvex constraints can be converted into linear constraints. In addition, we adopt the relaxation and penalty algorithm to deal with the Boolean constraint. Then, the optimization problem is equivalently trans-

#### **3. System Model and Problem Formulation** *3.1. System Model*

*Sensors* **2022**, *22*, x FOR PEER REVIEW 4 of 32

#### *3.1. System Model* As shown in Figure 1, we focus on the downlink of the massive MIMO LEO satellite

As shown in Figure 1, we focus on the downlink of the massive MIMO LEO satellite multigroup multicast communication system, and the LEO satellite uses a UPA, which is composed of *N* = *N<sup>x</sup>* × *N<sup>y</sup>* antennas and *L*(*L* ≤ *N*) RF links and covers *L* multicast groups and *K* active users, where *K* ≥ *L*, *K* ≥ *N*. Let the multicast group covered by the *l*th beam be *U<sup>l</sup>* and the number of users in this multicast group be |*U<sup>l</sup>* |, assuming that the number of users that can be accommodated in each DVB-S2X frame is *U<sup>s</sup>* . It should be noted that each user terminal belongs to only one multicast group, i.e., *U<sup>i</sup>* ∩ *U<sup>j</sup>* = ∅, ∀*i*, *j* ∈ {1, . . . , *L*}, *i* 6= *j*. Meanwhile, we assume that the user terminal is equipped with a single antenna capable of data stream demodulation. According to the DVB-S2X standard, in a transmission slot, multiple users' data in a multicast group are multiplexed into a specific forward error correction (FEC) codeword to provide services for more users. The service process is shown in Figure 2. multigroup multicast communication system, and the LEO satellite uses a UPA, which is composed of *N N N* = *x y* antennas and *L L N* ( ) RF links and covers *L* multicast groups and *K* active users, where *K L K N* , . Let the multicast group covered by the *l*th beam be *Ul* and the number of users in this multicast group be *Ul* , assuming that the number of users that can be accommodated in each DVB-S2X frame is *Us* . It should be noted that each user terminal belongs to only one multicast group, i.e., *U U i j L i j i j* = , , 1,...., , . Meanwhile, we assume that the user terminal is equipped with a single antenna capable of data stream demodulation. According to the DVB-S2X standard, in a transmission slot, multiple users' data in a multicast group are multiplexed into a specific forward error correction (FEC) codeword to provide services for more users. The service process is shown in Figure 2.

**Figure 1.** Transmission model of the massive MIMO LEO satellite multigroup multicast communication system. **Figure 1.** Transmission model of the massive MIMO LEO satellite multigroup multicast communication system.

**Figure 2.** Service process of the LEO communication satellite multigroup multicast transmission based on the DVB-S2X. **Figure 2.** Service process of the LEO communication satellite multigroup multicast transmission based on the DVB-S2X.

The received signal *kl*, *y* of the *k*th user in the *l*th multicast group can be expressed The received signal *yk*,*<sup>l</sup>* of the *k*th user in the *l*th multicast group can be expressed as:

as: , , , , [:, ] [:, ] , 1,..., , 1,..., *L H H k l k l RF BB l k j RF BB l k l <sup>l</sup> j j l <sup>y</sup> l s j s n k U l L* = <sup>=</sup> *h F F* <sup>+</sup> *h F F* + , (1) *yk*,*<sup>l</sup>* = *h H k*,*l FRFFBB*[:, *l*]*s<sup>l</sup>* + *L* ∑ *j*=1,*j*6=*l h H k*,*j FRFFBB*[:, *j*]*s<sup>l</sup>* + *nk*,*<sup>l</sup>* , *k* ∈ {1, . . . , |*U<sup>l</sup>* |}, *l* ∈ {1, . . . , *L*}, (1)

1, where the first term in (1) represents the expected received signal of the *k*th user in the *l*th multicast group, the second term represents the interference of other multicast groups and the third term represents the additive Gaussian white noise; 1 , *N k l h* represents the channel vector of the *k*th user in the *l*th multicast group, *L L BB F* represents the digital beamforming matrix, 1 [:, ] *<sup>L</sup> BB l F* represents the digital beamforming vector of the *l*th multicast group, and *N L RF F* represents the analog beamforming matrix, where each element of *FRF* should meet the unit modulus element [19], i.e., ( ) , 1 *RF i j F* = . In addition, *l s* represents the signal of the multicast group *Ul* , which meets the unit where the first term in (1) represents the expected received signal of the *k*th user in the *l*th multicast group, the second term represents the interference of other multicast groups and the third term represents the additive Gaussian white noise; *<sup>h</sup>k*,*<sup>l</sup>* <sup>∈</sup> <sup>C</sup>*N*×<sup>1</sup> represents the channel vector of the *<sup>k</sup>*th user in the *<sup>l</sup>*th multicast group, *<sup>F</sup>BB* <sup>∈</sup> <sup>C</sup>*L*×*<sup>L</sup>* represents the digital beamforming matrix, *<sup>F</sup>BB*[:, *<sup>l</sup>*] <sup>∈</sup> <sup>C</sup>*L*×<sup>1</sup> represents the digital beamforming vector of the *l*th multicast group, and *<sup>F</sup>RF* <sup>∈</sup> <sup>C</sup>*N*×*<sup>L</sup>* represents the analog beamforming matrix, where each element of *FRF* should meet the unit modulus element [19], i.e., (*FRF*)*i*,*<sup>j</sup>* <sup>=</sup> 1. In addition, *sl* represents the signal of the multicast group *U<sup>l</sup>* , which meets the unit power constraint, i.e., *E* n |*sl* | 2 o = 1, *n<sup>l</sup>* ∼ *CN*(0, *σ* 2 ) represents the additive Gaussian white noise, which is related to the Boltzmann constant *κ*, system bandwidth *B* and the noise temperature *T*.

power constraint, i.e., <sup>2</sup> 1 *<sup>l</sup>* = *s* , 2 (0, ) *<sup>l</sup> n CN* represents the additive Gaussian white noise, which is related to the Boltzmann constant , system bandwidth *B* and the noise temperature *T* . For the convenience of analysis, we set *<sup>F</sup>* <sup>∈</sup> <sup>C</sup>*N*×*<sup>L</sup>* <sup>=</sup> *<sup>F</sup>RFFBB* = [*<sup>f</sup>* 1 ,*f* 2 , . . . ,*f<sup>L</sup>* ] as the hybrid beamforming matrix, and *f <sup>l</sup>* <sup>∈</sup> <sup>C</sup>*N*×<sup>1</sup> <sup>=</sup> *<sup>F</sup>RFFBB*[:, *<sup>l</sup>*] is the hybrid beamforming vector of the *l*th multicast group. Therefore, (1) can be rewritten as:

$$y\_{k,l} = h\_{k,l}^H \mathbf{f}\_l \mathbf{s}\_l + \sum\_{j=1, j \neq l}^L h\_{k,j}^H \mathbf{f}\_j \mathbf{s}\_l + n\_{k,l}, k \in \{1, \dots, |\mathcal{U}\_l|\}, l \in \{1, \dots, L\} \tag{2}$$

, , , , 1, , 1,..., , 1,..., *L H H k l k l l l k j j l k l <sup>l</sup> j j l y s s n k U l L* = = + *h f h f* + , (2) Due to the high orbital speed of LEO satellites and the long transmission delay, it is Due to the high orbital speed of LEO satellites and the long transmission delay, it is difficult to obtain the precise instantaneous CSI. To cope with this problem, we adopt the statistical CSI, and the channel vector between the LEO satellite and the *k*th user in the *l*th multicast group at instant *t* and frequency *f* can be modeled as follows [20]:

difficult to obtain the precise instantaneous CSI. To cope with this problem, we adopt the

$$\mathbf{h}\_{k,l}(t,f) = \sum\_{p=1}^{P\_{k,l}} a\_{k,l,p} e^{j2\pi(f\_{d(k,l,p)}t - f\tau\_{k,l,p})} \times \mathbf{V}\_{k,l,p\_I} \tag{3}$$

( , , ) , , 2 ( ) , , , , , 1 ( , ) *k l d k l p k l p j f t f k l k l p k l p p t f a e* − = *h V* <sup>=</sup> , (3) where *f* denotes the carrier frequency, *k l p* , , *a* , *d k l p* ( , , ) *f* , *k l p* , , are the complex channel gain, where *f* denotes the carrier frequency, *ak*,*l*,*<sup>p</sup>* , *fd*(*k*,*l*,*p*) ,*τk*,*l*,*<sup>p</sup>* are the complex channel gain, Doppler shift and propagation delay, respectively, *Pk*,*<sup>l</sup>* denotes the number of propagation paths and *<sup>V</sup>k*,*l*,*<sup>p</sup>* <sup>∈</sup> <sup>C</sup>*N*×<sup>1</sup> is the UPA array response vector, which can be given by

,

*P*

$$\mathbf{V}\_{\mathbf{k},l,p} = \mathbf{V}(\boldsymbol{\varrho}\_{\mathbf{k},l,p}^{\mathbf{x}} \boldsymbol{\varrho}\_{\mathbf{k},l,p}^{\mathbf{y}}) = \boldsymbol{\sigma}\_{\mathbf{N}\_{\mathbf{x}}} \Big(\sin \boldsymbol{\varrho}\_{\mathbf{k},l,p}^{\mathbf{y}} \cos \boldsymbol{\varrho}\_{\mathbf{k},l,p}^{\mathbf{x}}\Big) \otimes \boldsymbol{\sigma}\_{\mathbf{N}\_{\mathbf{y}}} \Big(\cos \boldsymbol{\varrho}\_{\mathbf{k},l,p}^{\mathbf{y}}\Big),\tag{4}$$

*x y*

*k l p <sup>v</sup> <sup>v</sup>* , (4)

=

*l*th

$$\sigma\_{\rm Nr} \left( \sin \varphi\_{\rm k,l,p}^y \cos \varphi\_{\rm k,l,p}^x \right) = \frac{1}{\sqrt{N\_x}} \left( \mathbf{1}\_r e^{-j\frac{2\pi d}{\lambda} \sin \varphi\_{\rm k,l,p}^y \cos \varphi\_{\rm k,l,p}^x \dots} e^{-j\frac{2\pi d}{\lambda} (N\_\Gamma - 1) \sin \varphi\_{\rm k,l,p}^y \cos \varphi\_{\rm k,l,p}^x} \right) \tag{5}$$

$$
\sigma\_{N\_y}(\cos \varphi\_{k,l,p}^y) = \frac{1}{\sqrt{N\_y}} \left( 1, e^{-j\frac{2\pi d}{\lambda} \cos \varphi\_{k,l,p}^y}, \dots, e^{-j\frac{2\pi d}{\lambda} (N\_y - 1) \cos \varphi\_{k,l,p}^y} \right) \tag{6}
$$

where *ϕ x k*,*l*,*p* , *ϕ y k*,*l*,*p* represent azimuth angle and pitch angle associated with the propagation path *p* of the *k*th user in the *l*th multicast group, respectively, *λ* denotes the wavelength and *d* represents the spacing of antenna elements, the value of which is usually *λ*/2 [5].

Note that the LEO satellite communication system is usually operated under the line of sight (LOS) transmission, and the channel vector can be modeled using the widely accepted Rician distribution model as follows:

$$
\hbar\_{\mathbf{k},l} = \overline{\hbar}\_{\mathbf{k},l} + \overline{\hbar}\_{\mathbf{k},l} \tag{7}
$$

where *hk*,*<sup>l</sup>* = q*κk*,*lγk*,*<sup>l</sup> <sup>κ</sup>k*,*l*+<sup>1</sup> × *Vk*,*<sup>l</sup>* represents the LoS component, *h*e *<sup>k</sup>*,*<sup>l</sup>* = q *γk*,*<sup>l</sup> <sup>κ</sup>k*,*l*+<sup>1</sup> × *Vk*,*l*,*<sup>c</sup>* × *V H k*,*l* represents the multipath component, *<sup>κ</sup>k*,*<sup>l</sup>* denotes the Rician factor, *<sup>V</sup>k*,*l*,*<sup>c</sup>* <sup>∈</sup> <sup>C</sup>*Nut*×<sup>1</sup> <sup>∼</sup> CN(0, <sup>∑</sup>) represents Rician component, *Tr*(∑) = 1 and *γk*,*<sup>l</sup>* represents the average channel power, which mainly includes the transmit antenna gain *Gleo*, receiver antenna gain *Gut* and link power loss. The link power loss is mainly caused by the free space path loss *LPf s* and the atmospheric absorption loss *LPat*. Therefore, the average channel power *γk*,*<sup>l</sup>* can be written as

$$\gamma\_{k,l} = \mathbf{G}\_{l \text{elo}}[\mathbf{dB}] + \mathbf{G}\_{\text{ut}}[\mathbf{dB}] - LP\_{\text{at}}[\mathbf{dB}] - LP\_{fs}[\mathbf{dB}]\_{\text{t}} \tag{8}$$

where *LPf s* can be given by *LPf s* = 20 log10(*Dk*,*l*) + log10(*f*) + log10(4*π*/*c*), *c* is the speed of light, *Dk*,*<sup>l</sup>* represents the transmission distance, *LPat* is related to the carrier frequency, temperature *T*(*h*), pressure *P*(*h*) and humidity *ρ*(*h*), which can be given by *LPat* = R *<sup>h</sup>at hut LPat*(*f* , *T*(*h*), *P*(*h*), *ρ*(*h*))*dh*, *LPat*(*f* , *T*(*h*), *P*(*h*), *ρ*(*h*)) is the loss per meter, *hut* is the user's height and *hat* is the atmosphere thickness. The specific calculation method of *LPat* can be found in the literature [21].

For the convenience of analysis, it is assumed that the parameters in the channel vector *hk*,*<sup>l</sup>* are constant within coherence time and change over time in a certain ergodic process. In (3), the Doppler shift *fd*(*k*,*l*,*p*) and the propagation delay *τk*,*l*,*<sup>p</sup>* usually cause CSI errors. Next, we focus on analyzing the influence of propagation delay and Doppler shift on CSI errors.

*Doppler shift*: In the LEO satellite communication systems, the Doppler shift is usually large, which is mainly composed of the Doppler shift *f leo d*(*k*,*l*,*p*) generated by the LEO satellite motion and the Doppler shift *f ut d*(*k*,*l*,*p*) generated by the users' motion [22]. Since the transmission between the LEO satellite and user terminals is mainly under LOS, *f leo d*(*k*,*l*,*p*) of different transmission paths can be considered to be the same, and we omit the path index of *f leo d*(*k*,*l*,*p*) , i.e., <sup>n</sup> *f leo d*(*k*,*l*,*p*) o*P<sup>k</sup>* 1 = *f leo d*(*k*,*l*) ; *f leo d*(*k*,*l*) can be calculated using the LEO satellite ephemeris information and the location information of user terminals, as shown in Figure 3.

$$f\_{d(k,l)}^{\rm leo} = -\frac{f}{c} \times \frac{\omega^{\rm leo} r\_\varepsilon r \sin(\phi\_l - \phi\_{l\_0}) \mu(\theta\_{\rm max})}{\sqrt{r\_\varepsilon^2 + r^2 - 2r\_\varepsilon r \cos(\phi\_l - \phi\_{l\_0}) \mu(\theta\_{\rm max})}} \tag{9}$$

where *µ*(*θ*max) = cos- cos−<sup>1</sup> *re r* cos *θ*max <sup>−</sup> *<sup>θ</sup>*max , *r<sup>e</sup>* denotes the earth radius, *r* represents the distance between the LEO satellite track point and the earth center, (*φ<sup>t</sup>* − *φt*<sup>0</sup> ) represents the angular distance of the earth's surface along the LEO satellite trajectory from instant *t* to instant *t*0, *ω* represents the angular velocity of the LEO satellite and *c* represents the speed of light.

where

 

**Figure 3.** Geometric diagram of LEO satellite's orbital motion relative to the user terminal. **Figure 3.** Geometric diagram of LEO satellite's orbital motion relative to the user terminal.

( ) ( ) ( ) ( ) 0 0 max ( , ) 2 2 max sin 2 cos *leo e t t leo d k l <sup>e</sup> e t t r r <sup>f</sup> f c r r r r* − = − + − − , (9) The value of *f ut d*(*k*,*l*,*p*) in different transmission paths is different, which is mainly caused by the movement of user terminals and surrounding scatterers. The power spectrum of *f ut d*(*k*,*l*,*p*) follows the *Jakes* power spectrum model, and the normalized power spectrum can be expressed as:

$$S(f\_{d(k,l,p)}^{ut}) = \frac{1}{2\pi f\_{d(k,l,p),\max}^{ut}\sqrt{1 - \left(\frac{f\_{d(k,l,p)}^{ut}}{f\_{d(k,l,p),\max}^{ut}}\right)^2}}\tag{10}$$

the angular distance of the earth's surface along the LEO satellite trajectory from instant *t* to instant 0 *t* , represents the angular velocity of the LEO satellite and *c* represents the speed of light. The value of ( , , ) *ut d k l p f* in different transmission paths is different, which is mainly caused by the movement of user terminals and surrounding scatterers. The power spectrum of ( , , ) *ut d k l p f* follows the *Jakes* power spectrum model, and the normalized power spectrum can be expressed as: ( , , ) ( , , ) <sup>2</sup> ( , , ),max ( , , ),max 1 ( ) <sup>2</sup> 1 ( ) *ut d k l p ut ut d k l p d k l p ut d k l p S f f f f* = − , (10) The large Doppler shift in LEO satellite communication systems can make it difficult to receive correctly and result in the degradation of communication performance. In application, to mitigate the impact of Doppler shift, the solution of estimation and compensation is usually adopted [23]. The Doppler shift estimation mainly includes two steps: coarse estimation and fine estimation, which can refer to the literature [24]. When we obtain the estimated value of Doppler shift, *f e d* , the Doppler shift compensation of size *f cps <sup>d</sup>* = *f e d* can be implemented at the receivers. It should be noted that due to that the Doppler shift changes rapidly, in addition to the low SINR at the receivers and the limited pilot length in the DVB-S2X frame, the Doppler shift estimation is usually inaccurate, which can result in the incomplete compensation. Then, there would be the residual Doppler shift *f rsd d* , which can cause the sliding of channel phase. According to the Doppler shift estimation theory based the Cramer–Rao bound, the variance in the Doppler shift estimation can be expressed as the Cramer–Rao lower bound (CRLB) [25], i.e.,

$$
\sigma\_{f\_d^\varepsilon}^2 = \text{CRLB}(f) = \frac{1}{\text{SNR}} \frac{3}{2\pi^2 T^2 N (N^2 - 1)}\tag{11}
$$

sation is usually adopted [23]. The Doppler shift estimation mainly includes two steps: coarse estimation and fine estimation, which can refer to the literature [24]. When we obtain the estimated value of Doppler shift, *e d f* , the Doppler shift compensation of size *cps e d d f f* = can be implemented at the receivers. It should be noted that due to that the Dopwhere *N* is the pilot length, *T* is the sampling time and SNR is the signal-to-noise ratio. According to the properties of variance, after Doppler shift compensation of size *f cps d* , the variance of residual Doppler shift *f rsd d* is equal to the variance of *f e d* , i.e.,

$$
\sigma\_{f\_d^{\text{rad}}}^2 = \sigma\_{f\_d^e}^2 = \frac{1}{\text{SNR}} \frac{3}{2\pi^2 T^2 N (N^2 - 1)}\tag{12}
$$

which can result in the incomplete compensation. Then, there would be the residual Doppler shift *rsd d f* , which can cause the sliding of channel phase. According to the Doppler shift estimation theory based the Cramer–Rao bound, the variance in the Doppler shift estimation can be expressed as the Cramer–Rao lower bound (CRLB) [25], i.e., The influence of residual Doppler shift on channel phase errors can be expressed as *φf rsd d* = 2*π f rsd d* ∆*T*, where ∆*T* represents the sum of the downlink propagation delay and the duration of the DVB-S2X frame. Therefore, the variance of *φ<sup>f</sup> rsd d* can be expressed as *σ* 2 *φ f rsd d* = 4*π* 2*σ* 2 *f rsd d* ∆*T* <sup>2</sup> = <sup>1</sup> *SNR* 6∆*T* 2 *<sup>T</sup>*2*N*(*N*2−1) , and *φ<sup>f</sup> rsd d* follows a real-valued Gaussian distribution with mean zero and variance *σ* 2 *φ f rsd d* , i.e., *φ<sup>f</sup> rsd d* ∼ *N* 0, *σ* 2 *φ f rsd d* . The channel error vector caused by *φ<sup>f</sup> rsd* can be expressed as *v<sup>φ</sup>* = [*e jφ f rsd d* ,1 ,*e jφ f rsd d* ,2 , . . . ,*e jφ f rsd d* ,*N* ] *T* .

*rsd*

*d f d Propagation Delay*: The orbital height of LEO satellites is about 300 km to 2000 km, and the long transmission distance can cause the larger propagation delay. Note that the influence of atmosphere on propagation delay mainly includes ionospheric delay and tropospheric delay [26]. When the signal passes through the ionosphere, due to the refraction effect of the electromagnetic wave, the propagation path and speed of the signal will change. Meanwhile, the ionospheric delay is irregular, which is difficult to describe with a physical model. When the signal passes through the troposphere, the propagation speed, direction and path of the signal will change, which can result in propagation delay. The tropospheric delay is related to air pressure, air humidity and satellite elevation. The commonly used tropospheric delay correction model is given in the literature [27,28]. The round-trip delay of the LEO satellite with an orbit altitude of 1200 km is about 20 ms. The long propagation delay will lead to the expiration of CSI, which can result in CSI errors, phase disturbance and other problems. To handle this problem, the delay compensation of size *τ cps* <sup>=</sup> *βτ*min + (<sup>1</sup> <sup>−</sup> *<sup>β</sup>*)*τ*max, (<sup>0</sup> <sup>≤</sup> *<sup>β</sup>* <sup>≤</sup> <sup>1</sup>) is usually implemented at the receivers, and *τ* min *<sup>k</sup>*,*<sup>l</sup>* <sup>=</sup> min<sup>n</sup> *τk*,*l*,*<sup>p</sup>* o*Pk*,*<sup>l</sup>* 1 and *τ* max *<sup>k</sup>*,*<sup>l</sup>* <sup>=</sup> max<sup>n</sup> *τk*,*l*,*<sup>p</sup>* o*Pk*,*<sup>l</sup>* 1 represent the minimum propagation delay and the maximum propagation delay of the *k*th user in the *l*th multicast group, respectively.

However, due to that the atmospheric propagation delay is irregular, the transmission delay cannot be fully compensated. Therefore, the incomplete delay compensation, expired CSI and distortion of high-frequency devices would cause the channel phase disturbance [29]. Let the phase disturbance be *φτ*, which follows a real-valued Gaussian distribution with mean zero and variance *σ* 2 *φτ* , i.e., *φ<sup>τ</sup>* ∼ *N* 0, *σ* 2 *φτ* . The channel error vector caused by *φ<sup>τ</sup>* can be expressed as *vφ<sup>τ</sup>* = [*e <sup>j</sup>φτ*,1 ,*e <sup>j</sup>φτ*<sup>2</sup> , . . . ,*e <sup>j</sup>φτ*,*<sup>N</sup>* ] *T* .

In conclusion, considering the influence of the residual Doppler shift and the phase disturbance, the relationship between the real channel vector *hk*,*<sup>l</sup>* and the estimated channel vector *h*ˆ *<sup>k</sup>*,*<sup>l</sup>* can be expressed as:

$$\mathfrak{h}\_{k,l} = \mathfrak{h}\_{k,l} \odot \mathfrak{v}\_{\Phi\_{\tilde{d},k,l}} \odot \mathfrak{v}\_{\Phi\_{\tilde{r},k,l}} = \operatorname{diag}\left(\operatorname{diag}\left(\mathfrak{h}\_{k,l}\right) \mathfrak{v}\_{\Phi\_{\tilde{d},k,l}}\right) \mathbf{v}\_{\Phi\_{\tilde{k},l'}} \tag{13}$$

where  represents the Hadamard product. Let the channel phase of the *k*th user in the *l*th multicast group be *θk*,*<sup>l</sup>* = [*θk*,*l*,1, *θk*,*l*,2, . . . , *θk*,*l*,*N*] *T* , which satisfies the uniform distribution between 0 ∼ 2*π*. Then, the real channel phase with phase errors at instant *t*<sup>1</sup> is as follows:

$$
\theta\_{k,l}(t\_1) = \theta\_{k,l}(t\_0) + \phi\_{f\_{d,k,l}^{\rm rad}} + \phi\_{\sf T\_{k,l}} \tag{14}
$$

#### *3.2. Problem Formulation*

#### 3.2.1. User Clustering

Before the joint user scheduling and hybrid beamforming design, it is necessary to group the active users within the coverage of the LEO satellite. Based on the CSI, the hierarchical clustering algorithm is adopted to group users [30]. As shown in Figures 4 and 5, the hierarchical clustering algorithm adopts the bottom-up method, where each user initially forms a group, and then according to the similarity measurement function, the user groups which meet the similarity threshold constraint are combined until the desired number of groups is formed.

the *l*th

follows:

the *l*th

follows:

tion between

tion between

*3.2. Problem Formulation* 3.2.1. User Clustering

*3.2. Problem Formulation* 3.2.1. User Clustering

number of groups is formed.

number of groups is formed.

multicast group be

0 2

multicast group be

0 2

where represents the Hadamard product. Let the channel phase of the

 

*k l k l k l k l N* =

, , ,1 , ,2 , , [ , ,..., ]*<sup>T</sup>*

 

 

where represents the Hadamard product. Let the channel phase of the

*Sensors* **2022**, *22*, x FOR PEER REVIEW 9 of 32

, , ,1 , ,2 , , [ , ,..., ]*<sup>T</sup>*

 

 

, 1 , 0 ( ) ( ) *rsd k l d k l k l k l <sup>f</sup> t t*

 

= + +

, 1 , 0 ( ) ( ) *rsd k l d k l k l k l <sup>f</sup> t t*

Before the joint user scheduling and hybrid beamforming design, it is necessary to group the active users within the coverage of the LEO satellite. Based on the CSI, the hierarchical clustering algorithm is adopted to group users [30]. As shown in Figures 4 and 5, the hierarchical clustering algorithm adopts the bottom-up method, where each user initially forms a group, and then according to the similarity measurement function, the user groups which meet the similarity threshold constraint are combined until the desired

Before the joint user scheduling and hybrid beamforming design, it is necessary to group the active users within the coverage of the LEO satellite. Based on the CSI, the hierarchical clustering algorithm is adopted to group users [30]. As shown in Figures 4 and 5, the hierarchical clustering algorithm adopts the bottom-up method, where each user initially forms a group, and then according to the similarity measurement function, the user groups which meet the similarity threshold constraint are combined until the desired

= + +

 

. Then, the real channel phase with phase errors at instant

. Then, the real channel phase with phase errors at instant

 

, ,

, ,

,

,

 

,

,

 

 

*k l k l k l k l N* =

 

*k*th

1 *t*

*k*th

, which satisfies the uniform distribu-

, which satisfies the uniform distribu-

user in

is as

(14)

is as

(14)

1 *t*

user in

**Figure 4.** Schematic diagram of the hierarchical clustering algorithm. **Figure 4.** Schematic diagram of the hierarchical clustering algorithm.

**Figure 4.** Schematic diagram of the hierarchical clustering algorithm.

**Figure 5. Figure 5.** Schematic diagram of multibeam coverage with 7 multicast groups. Schematic diagram of multibeam coverage with 7 multicast groups.

**Figure 5.** Schematic diagram of multibeam coverage with 7 multicast groups. We adopt the similarity measurement function among multicast groups based on the *Ward* connection method, i.e.,

$$d(i,j) = \sqrt{\frac{2n\_in\_j}{n\_i + n\_j}} dist(h\_i^{eq}, h\_j^{eq}),\tag{15}$$

where *n<sup>i</sup>* , *n<sup>j</sup>* represent the number of users of the group *i* and the group *j*, respectively, *h eq i* , *h eq j* represent the equivalent CSI of the group *i* and the group *j*, respectively, and *dist*(*h eq i* , *h eq j* ) represents the Euclidean distance between the vector *h eq i* and the vector *h eq i* , i.e.,

$$dist(\mathsf{h}\_i^{eq}, \mathsf{h}\_j^{eq}) = \left\| \frac{\mathsf{h}\_i^{eq}}{\left\| \mathsf{h}\_i^{eq} \right\|} - \frac{\mathsf{h}\_j^{eq}}{\left\| \mathsf{h}\_j^{eq} \right\|} \right\|\tag{16}$$

3.2.2. System Rate

Affected by CSI errors, both the ergodic communication rate and the ergodic SINR do not admit explicit expressions. To handle this challenge, the statistical average method is adopted to model the SINR and the communication rate. Therefore, the SINR and the communication rate of the *k*th user in the *l*th multicast group can be expressed as:

$$\mathbb{E}\{\text{SINR}\_{k,l}\} \approx \frac{\mathbb{E}\left\{ \left| h\_{k,l}^H f\_l \right|^2 \right\}}{\mathbb{E}\left\{ \sum\_{j=1, j\neq l}^L \left| h\_{k,l}^H f\_j \right|^2 \right\} + \sigma^2},\tag{17}$$

$$R\_{k,l} \approx B \log\_2 \left( 1 + \frac{\mathbb{E}\left\{ \left| h\_{k,l}^H \boldsymbol{\mathcal{f}}\_l \right|^2 \right\}}{\mathbb{E}\left\{ \sum\_{j=1, j \neq l}^L \left| h\_{k,l}^H \boldsymbol{\mathcal{f}}\_j \right|^2 \right\} + \sigma^2} \right), \tag{18}$$

where *B* denotes system bandwidth. Equations (17) and (18) are approximations with closed form, the feasibility of which have been discussed in detail in Refs. [31,32].

#### 3.2.3. Problem Description

We take the system EE as the optimization objective, and the EE is defined as the ratio of the system communication rate to the total power consumption, which can be modeled as:

$$\text{EE} = \frac{\sum\_{l=1}^{L} \min\left( \left\{ R\_{k,l} \right\}\_{k=1}^{\left| U\_l \right|} \right)}{P\_{\text{total}}} = \frac{B \sum\_{l=1}^{L} \log\_2 \left( 1 + \min\left\{ \frac{\left| h\_k^H f\_l \right|^2}{\sum\_{j=1, j \neq l}^{L} \left| h\_k^H f\_j \right|^2 + \sigma^2} \right\}\_{k=1}^{\left| U\_s \right|} \right)}{P\_t + P\_0},\tag{19}$$

where *Ptotal* represents the total power consumption, *P<sup>t</sup>* = *L* ∑ *l f l f H l*  denotes the transmission power of the LEO satellite, and *P*<sup>0</sup> denotes the inherent power consumption of the communication system.

Let the Boolean variable *ηk*,*<sup>l</sup>* ∈ {0, 1} indicate whether the *k*th user in the *l*th multicast group is served, *ηk*,*<sup>l</sup>* = 1 and *ηk*,*<sup>l</sup>* = 0 indicate that the user can be served and not served, respectively, and *η* = [*η*1, *η*2, . . . , *ηL*], *η<sup>l</sup>* = h *η*1,*<sup>l</sup>* , *η*2,*<sup>l</sup>* , . . . , *<sup>η</sup>*|*U<sup>l</sup>* |,*L* i*T* . In conclusion, under the constraints of the transmission power and QoS, the problem of maximizing system EE can be modeled as:

$$Q\_1: \max\_{\eta f\_l \in \text{SINR}\_l^{\text{min}}} \text{EE} = \frac{B \sum\_{l=1}^L \log\_2 \left( 1 + \text{SINR}\_l^{\text{min}} \right)}{\sum\_{l}^L \left| f\_l f\_l^H \right|^2 + P\_0},\tag{20}$$

$$s.t. \ C\_1: \eta\_{k,l} \in \{0, 1\}, \forall k, l,\tag{21}$$

$$\mathbf{C}\_2: \mathbf{SINR}\_{k,l} \ge \eta\_{k,l} \mathbf{SINR}\_l^{\min}, \forall k, l\_\prime \tag{22}$$

$$\mathcal{C}\_3: \text{SINR}\_l^{\text{min}} \ge \text{SINR}\_{0\prime} \,\forall l \,\, \tag{23}$$

$$\mathbb{C}\_{4}: \sum\_{k=1}^{\binom{|\mathcal{U}\_{l}|}{2}} \eta\_{k,l} = \mathbb{U}\_{\mathbf{s}\_{\mathsf{V}}} \forall l\_{\mathsf{\cdot}} \tag{24}$$

$$\mathbf{C}\_{\mathsf{S}} : \sum\_{l=1}^{L} \left| f\_{l} f\_{l}^{H} \right|^{2} \leq P\_{\mathsf{T}\_{\mathsf{V}}} \tag{25}$$

where SINRmin *l* represents the minimum SINR of the *l*th multicast group and constraint *C*<sup>3</sup> represents that SINRmin *l* should be greater than the minimum SINR constraint SINR0. Constraint *C*<sup>4</sup> limits the number of scheduled users in each multicast group to *U<sup>s</sup>* .

#### **4. Joint User Scheduling and Hybrid Beamforming Design for Maximizing EE**

In this section, we focus on the robust joint user scheduling and hybrid beamforming design strategy to maximize system EE. To handle the QCQP form problem and nonconvexity in optimization problem *Q*1, the SDP method is applied to make the optimization problem more tractable. Then, we transform the optimization problem *Q*<sup>1</sup> into a DC programming problem. To address the DC programming problem, we adopt the CCCP algorithm. Finally, a penalty iterative algorithm is adopted to handle the rank-one matrix constraint.

#### *4.1. SDP Algorithm*

It is worth noting that the objective function and constraints *C*<sup>2</sup> and *C*<sup>5</sup> in the problem *Q*<sup>1</sup> involve the quadratic form of the variable *f l* , therefore, *Q*<sup>1</sup> is the QCQP form problem. To handle this problem, we invoke the SDP algorithm, a new variable *W<sup>l</sup>* , *f l f H l* is introduced, and the positive semidefinite matrix *W<sup>l</sup>* needs to meet the constraints of *W<sup>l</sup>*0 and *rank*(*Wl*) = 1. Then, the problem *Q*<sup>1</sup> can be equivalent to:

$$Q\_2: \max\_{\eta, \mathbf{W}\_l, \text{SINR}\_l^{\text{min}}} \text{EE} = \frac{B \sum\_{l=1}^L \log\_2 \left( 1 + \text{SINR}\_l^{\text{min}} \right)}{\sum\_{l}^L \text{Tr}(\mathbf{W}\_l) + P\_0} \tag{26}$$

$$\text{s.t.}\ \mathsf{C}\_{1}\mathsf{C}\_{2}\ \mathsf{C}\_{3}\ \mathsf{C}\_{4}\ \mathsf{in}\ \mathsf{Q}\_{1}\!\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/\!/]}$$

$$\mathbb{C}\_{5}: \sum\_{l}^{L} \text{Tr}(\mathbf{W}\_{l}) \le \mathbf{P}\_{\mathsf{T}\_{l}} \tag{28}$$

$$\mathsf{C}\_{6}: \mathsf{W}\_{l} \succeq \mathsf{0}, \;\forall l,\tag{29}$$

$$\mathbf{C}\_{\mathcal{T}} : \operatorname{rank}(\mathbf{W}\_l) = 1, \,\forall l,\,\tag{30}$$

Similarly, the SINR and the communication rate of the *k*th user in the *l*th multicast group can be equivalently converted to:

$$\text{SINR}\_{k,l} \approx \frac{\text{E}\left\{ \text{Tr}(\mathbf{H}\_{k,l}\mathbf{W}\_l) \right\}}{\text{E}\left\{ \sum\_{j=1, j\neq l}^{L} \text{Tr}\left(\mathbf{H}\_{k,l}\mathbf{W}\_j\right) \right\} + \sigma^2},\tag{31}$$

$$\begin{split} R\_{\mathbb{k}J} &= B \log\_2 \left( 1 + \frac{\mathbb{E} \left\{ \text{Tr} \left( \mathbf{H}\_{\mathbb{k}J} \mathbf{W}\_l \right) \right\}}{\mathbb{E} \left\{ \sum\_{j=1, j \neq l}^{\mathbb{L}} \text{Tr} \left( \mathbf{H}\_{\mathbb{k}J} \mathbf{W}\_j \right) \right\} + \sigma^2} \right) = B \log\_2 \left( 1 + \frac{\text{Tr} \left\{ \mathbb{E} \left\{ \mathbf{H}\_{\mathbb{k}J} \mathbf{W}\_l \right\} \right\}}{\sum\_{j=1, j \neq l}^{\mathbb{L}} \text{Tr} \left\{ \mathbb{E} \left\{ \mathbf{H}\_{\mathbb{k}J} \mathbf{W}\_l \right\} \right\} + \sigma^2} \right) \\ &= B \log\_2 \left( 1 + \frac{\text{Tr} \left( \mathbf{H}\_{\mathbb{k}J} \mathbf{W}\_l \right)}{\sum\_{j=1, j \neq l}^{\mathbb{L}} \text{Tr} \left( \mathbf{H}\_{\mathbb{k}J} \mathbf{W}\_j \right) + \sigma^2} \right) \end{split} \tag{32}$$

where *Hk*,*<sup>l</sup>* ∈ *C <sup>M</sup>*×*<sup>M</sup>* is the instantaneous channel autocorrelation matrix of the *k*th user in the *l*th multicast group and *Hk*,*<sup>l</sup>* ∈ *C <sup>M</sup>*×*<sup>M</sup>* is the long-term channel autocorrelation matrix. The relationship between the two can be expressed as:

$$\overline{\mathbf{H}}\_{k,l} = \mathbb{E}\left\{ \mathbf{H}\_{k,l} \right\} \stackrel{\Delta}{=} \mathbb{E}\left\{ \mathbf{\hat{h}}\_{k,l} \mathbf{\hat{h}}\_{k,l}^H \right\} = \text{diag}\left( \mathbf{\hat{h}}\_{k,l} \right) \mathbf{P}\_{f\_{d,k}^{\text{red}}} \mathbf{Q}\_{\mathbf{\tau}\_{k,l}} \text{diag}\left( \mathbf{\hat{h}}\_{k,l}^H \right), \tag{33}$$

where *P<sup>f</sup> rsd d*,*k*,*l* and *Qτk*,*<sup>l</sup>* can be expressed as follows:

*Pf rsd d*,*k*,*l* = *E* ( *vφ f rsd d*,*k*,*l v H φ f rsd d*,*k*,*l* ) = E ( *e jφ f rsd d*,*k*,*l* ,1 ,*e jφ f rsd d*,*k*,*l* ,2 , . . . ,*e jφ f rsd d*,*k*.*l* ,*M T e* −*jφ f rsd d*,*k*,*l* ,1 ,*e* −*jφ f rsd d*,*k*,*l* ,2 , . . . ,*e* −*jφ f rsd d*,*k*.*l* ,*M* ) = E 1 · · · *e jφ f rsd d*,*k*,*l* ,1 *e* −*jφ f rsd d*,*k*.*l* ,*M* . . . . . . . . . *e jφ f rsd d*,*k*.*l* ,*M e* −*jφ f rsd d*,*k*,*l* ,1 · · · 1 = 1 · · · *E e jφ f rsd d*,*k*,*l* ,1 *e* −*jφ f rsd d*,*k*.*l* ,*M* . . . . . . . . . *E e jφ f rsd d*,*k*.*l* ,*M e* −*jφ f rsd d*,*k*,*l* ,1 · · · 1 (34)

In (34), the diagonal elements of *P<sup>f</sup> rsd d*,*k*,*l* are all 1, and the elements in row *i* and column *j* on the non-diagonal are E *e jφ f rsd d*,*k*.*l* ,*i e* −*jφ f rsd d*,*k*,*l* ,*j* = E *e jφ f rsd d*,*k*.*l* ,*i* E *e* −*jφ f rsd d*,*k*,*l* ,*j* , according to *φf rsd d* ∼ *N* 0, *σ* 2 *φ f rsd d* , *E jφ f rsd* ,*i jφ f rsd* 1 − *φ* 2 *f rsd d*,*k*,*l* 2*σ* 2 *φ f rsd*

$$\begin{split} \left(e^{\frac{\beta\theta\_{f}f\_{d,l}^{sd}}{\sqrt{d}}}\right)^{} &= \int\_{-\infty}^{\infty} e^{\frac{\beta\theta\_{f}f\_{d,l}^{sd}}{\sqrt{2\pi}\sigma\_{\theta\_{f}f\_{d,l}^{sd}}}} \frac{1}{\sqrt{2\pi}\sigma\_{\theta\_{f}f\_{d,l}^{sd}}} e^{-\frac{\beta\theta\_{f}f\_{d}^{sd}}{\sqrt{d}}} d\Phi\_{f\_{d,l}^{sd},i} \\ &= e^{\frac{\beta\theta\_{f}f\_{d,l}^{sd}}{2}} \int\_{-\infty}^{\infty} \frac{1}{\sqrt{2\pi}\sigma\_{\theta\_{f}f\_{d,l}^{sd}}} e^{-\frac{\beta\theta\_{f}f\_{d,l}^{sd}}{2\sigma\_{\theta\_{f}f\_{d,l}^{sd}}^{2}}} d\Phi\_{f\_{d,l}^{sd},i} \\ &= e^{2\beta} \end{split} \tag{35}$$

$$\begin{split} \text{Similarly, } & \mathbb{E}\left\{ e^{-j\Phi\_{f^{cd}\_{d\lambda},j}} \right\} = e^{-\frac{\sigma\_{\tilde{f}^{cd}\_{d}}}{2}}. \text{ Therefore, } \mathbb{E}\left\{ e^{j\Phi\_{f^{cd}\_{d\lambda},j}} e^{-j\phi\_{f^{cd}\_{d\lambda},j}} \right\} = \mathbb{E}\left\{ e^{j\frac{\phi\_{f^{cd}\_{d\lambda},j}}{2}} \right\}, \\ & \mathbb{E}\left\{ e^{-j\phi\_{f^{cd}\_{d\lambda},j}} \right\} = e^{-\sigma\_{\tilde{f}^{cd}\_{d\lambda}}}. \end{split}$$

$$\begin{aligned} \mathbf{Q}\_{\mathbb{R},j} &= \mathrm{E}\left\{ \mathbf{z}\_{\theta\_{\mathbb{R},j}} \boldsymbol{\sigma}\_{\theta\_{\mathbb{R},j}}^{H} \right\} \\ &= \mathrm{E}\left\{ \left[ e^{j\Phi\_{\mathbb{R},j}}, e^{j\Phi\_{\mathbb{R},j}}, \dots, e^{j\Phi\_{\mathbb{R},j},M} \right]^{T} \left[ e^{-j\Phi\_{\mathbb{R},j},1}, e^{-j\Phi\_{\mathbb{R},j},2}, \dots, e^{-j\Phi\_{\mathbb{R},j},M} \right] \right\} \\ &= \mathrm{E}\left\{ \begin{array}{cccc} 1 & \cdots & \cdots & e^{j\Phi\_{\mathbb{R},j},1}e^{-j\Phi\_{\mathbb{R},j},M} \\ \vdots & \ddots & \vdots \\ e^{j\Phi\_{\mathbb{R},j},M}e^{-j\Phi\_{\mathbb{R},j},1} & \cdots & 1 \\ \end{array} \right\} \\ &= \left\{ \begin{array}{cccc} 1 & \cdots & \cdots & E\left\{ e^{j\Phi\_{\mathbb{R},j},1}e^{-j\Phi\_{\mathbb{R},j},M} \right\} \\ \vdots & \ddots & \vdots \\ E\left\{ e^{j\Phi\_{\mathbb{R},j},M}e^{-j\Phi\_{\mathbb{R},j},1} \right\} & \cdots & 1 \end{array} \right\} \end{aligned} \tag{36}$$

In (36), the diagonal elements of *Qτk*,*<sup>l</sup>* are all 1, and the elements in row *i* and column *j* on the non-diagonal are E n *e <sup>j</sup>φτk*,*<sup>l</sup>* ,*i e* −*jφτk*,*<sup>l</sup>* ,*j* o = E n *e <sup>j</sup>φτk*,*<sup>l</sup>* ,*i* o E n *e* −*jφτk*,*<sup>l</sup>* ,*j* o , according to *φ<sup>τ</sup>* ∼ *N* 0, *σ* 2 *φτ* ,

$$\begin{split} \mathrm{E}\left\{ e^{j\phi\_{\mathbb{T}\_{k,l},i}} \right\} &= \int\_{-\infty}^{\infty} e^{j\phi\_{\mathbb{T}\_{k,l},i}} \frac{1}{\sqrt{2\pi}\sigma\_{\phi\_{\mathbb{T}\_{k,l}}}} e^{-\frac{\phi\_{\mathbb{T}\_{k,l},i}^{2}}{2\sigma\_{\phi\_{\mathbb{T}\_{k,l}}}^{2}}} d\phi\_{\mathbb{T}\_{k,l},i} \\ &= e^{-\frac{\sigma\_{\phi\_{\mathbb{T}\_{k,l}}}^{2}}{2}} \int\_{-\infty}^{\infty} \frac{1}{\sqrt{2\pi}\sigma\_{\phi\_{\mathbb{T}\_{k,l}}}} e^{-\frac{\left(\Phi\_{\mathbb{T}\_{k,l},i} - j\sigma\_{\phi\_{\mathbb{T}\_{k,l}}}\right)^{2}}{2\sigma\_{\phi\_{\mathbb{T}\_{k,l}}}^{2}}} d\phi\_{\mathbb{T}\_{k,l},i} \\ &= e^{-\frac{\sigma\_{\phi\_{\mathbb{T}\_{k,l}}^{2}}^{2}}} \end{split} \tag{37}$$

$$\begin{split} & \text{Similarly,} \operatorname{\mathbf{E}} \left\{ e^{j\phi\_{\tau\_{\mathbf{k},l}}} \right\} = e^{-\frac{\sigma\_{\theta\tau\_{\mathbf{k},l}}^{2}}{2}}. \operatorname{\mathbf{Therefore,}} \operatorname{\mathbf{E}} \left\{ e^{j\phi\_{\tau\_{\mathbf{k},l}}} e^{-j\phi\_{\tau\_{\mathbf{k},l},i}} \right\} = \operatorname{\mathbf{E}} \left\{ e^{j\phi\_{\tau\_{\mathbf{k},l}}} \right\} \operatorname{\mathbf{E}} \left\{ e^{-j\phi\_{\tau\_{\mathbf{k},l},i}} \right\} \\ & = e^{-\sigma\_{\phi\_{\tau\_{\mathbf{k},l}}^{2}}}. \end{split}$$

#### *4.2. DC Programming*

Since the constraint *C*<sup>1</sup> is a Boolean constraint and the constraint *C*<sup>2</sup> is a nonconvex constraint, the problem *Q*<sup>2</sup> is a nonconvex and nonsmooth combinatorial optimization problem. To handle this challenge, we can transform the problem *Q*<sup>2</sup> into a DC programming problem [33]. Therefore, the relaxation variable *ζk*,*<sup>l</sup>* is introduced as the lower bound of the SINR of the *<sup>k</sup>*th user in the *<sup>l</sup>*th multicast group, *<sup>ζ</sup>* = [*ζ*1, *<sup>ζ</sup>*2, . . . , *<sup>ζ</sup>L*], *<sup>ζ</sup><sup>l</sup>* = [*ζl*,1, *<sup>ζ</sup>l*,2, . . . , *<sup>ζ</sup>l*,|*U<sup>l</sup>* | ] *T* . The problem *Q*<sup>2</sup> can be equivalently converted to:

$$Q\_3: \max\_{\eta, \mathbf{W}, \text{SINR}^{\text{min}}\_l, \mathcal{J}} \text{EE} = \frac{B \sum\_{l=1}^L \log\_2 \left( 1 + \text{SINR}^{\text{min}}\_l \right)}{\sum\_{l=1}^L \text{Tr}(\mathbf{W}\_l) + P\_0},\tag{38}$$

$$\text{s.t.}\,\mathsf{C}\_{1}\,\mathsf{C}\_{3}\,\mathsf{C}\_{4}\,\mathsf{C}\_{5}\,\mathsf{C}\_{6}\,\mathsf{C}\_{6}\,\mathsf{C}\_{7}\,\text{in}\,\mathsf{Q}\_{1},\tag{39}$$

$$\mathsf{C}\_{2}: \mathsf{SINR}\_{k,l} \geq \mathcal{J}\_{k,l} \,\forall k, l,\tag{40}$$

$$\mathsf{C}\_{8}: \mathbb{zeta}\_{k,l} \geq \eta\_{k,l} \mathbf{SINR}\_{l}^{\min}, \forall k, l,\tag{41}$$

where the constraint *C*<sup>2</sup> can be equivalently converted to:

$$\mathbf{C}\_2 \Rightarrow \mathbf{1} + \mathbf{SINR}\_{k,l} \ge \mathbf{1} + \mathbf{\zeta}\_{k,l\prime} \tag{42}$$

To further express (42) in the form of DC programming, we introduce new function variables Γ*k*,*l*(*W*) and I*k*,*l*(*W*, *ζk*,*l*):

$$\Gamma\_{k,l}(\mathbf{W}) = \sigma^2 + \sum\_{j=1, j \neq l}^{L} \text{Tr}\left(\mathbf{H}\_{k,l}\mathbf{W}\_l\right), \tag{43}$$

$$\mathbf{I}\_{k,l}(\mathbf{W}, \tilde{\xi}\_{k,l}) = \frac{\sigma^2 + \sum\_{j=1}^{L} \text{Tr}\left(\mathbf{H}\_{k,l} \mathbf{W}\_j\right)}{1 + \tilde{\xi}\_{k,l}},\tag{44}$$

Therefore, the constraint *C*<sup>2</sup> can be rewritten as:

$$\mathbf{C}\_{2} \Rightarrow \Gamma(\mathbf{W}) - \mathbf{I}(\mathbf{W}\_{\prime}\mathbb{Z}\_{k,l}) \leq \mathbf{0},\tag{45}$$

In (45), Γ*k*,*l*(*Wl*) is the affine function of *W*, *Ik*,*l*(*W*, *ζk*,*l*) is the concave function of *W* and *ζk*,*<sup>l</sup>* . The transformed constraint *C*<sup>2</sup> is a typical DC constraint.

Similarly, the constraint *C*<sup>8</sup> can be equivalently converted into the following DC form:

$$\mathbf{C}\_{8} \Rightarrow \mathbf{4}\zeta\_{k,l} + \left(\eta\_{k,l} - \text{SINR}\_{l}^{\text{min}}\right)^{2} \geq \left(\eta\_{k,l} + \text{SINR}\_{l}^{\text{min}}\right)^{2},\tag{46}$$

In (38), the objective function in the problem *Q*<sup>3</sup> is a fractional programming problem with the sum-of-ratios form. To handle this problem, we invoke the quadratic transformation algorithm [34] and convert the problem *Q*<sup>3</sup> into the following form:

$$Q\_4: \max\_{\eta, \mathcal{W}, \text{SINR}^{\text{min}}\_l} \text{EE} = 2q \left( B \sum\_{l=1}^L \mathcal{Y}\_l(\text{SINR}^{\text{min}}\_l) \right)^{\frac{1}{2}} - q^2 \left( \sum\_{l=1}^L \text{Tr}(\mathcal{W}\_l) + P\_0 \right), \tag{47}$$

*s*.*t*. *C*1, *C*2, *C*3, *C*4, *C*<sup>5</sup> ,*C*<sup>6</sup> , *C*7, *C*<sup>8</sup> *in Q*3, (48)

where *q* is the introduced auxiliary variable, and *Υ<sup>l</sup>* (SINRmin *l* ) is the introduced auxiliary function, which can be expressed as:

$$Y\_I(\text{SINR}\_l^{\text{min}}) = \log\_2\left(1 + \text{SINR}\_l^{\text{min}}\right),\tag{49}$$

$$q = \frac{\sqrt{\sum\_{l=1}^{L} \mathcal{Y}\_l(\text{SINR}\_l^{\text{min}})}}{\sum\_{l=1}^{L} \text{Tr}(\mathbf{W}l) + P\_0} \text{ }\tag{50}$$

In addition, for the nonsmooth combinatorial optimization problem caused by constraint *C*1, we invoke a relaxation and penalty algorithm. Firstly, we relax constraint *C*<sup>1</sup> into *C*<sup>1</sup> ⇒ 0 ≤ *ηk*,*<sup>l</sup>* ≤ 1, ∀*k*, *l* . Meanwhile, to avoid the non-duality of the solution of *ηk*,*<sup>l</sup>* caused by the relaxation, the penalty term P(*ηk*,*l*) = *ηk*,*<sup>l</sup>* log *ηk*,*<sup>l</sup>* + (1 − *ηk*,*l*)log(1 − *ηk*,*l*) is introduced into the objective function: let *λ*<sup>1</sup> > 0 be the penalty factor, and the problem *Q*<sup>4</sup> can be equivalently converted to:

$$Q\_5: \max\_{\eta, \mathsf{W} \text{SINR}^{\text{min}}\_l, \mathcal{J}} \mathsf{EE} = 2q \left( B \sum\_{l=1}^L Y\_l(\mathsf{SINR}^{\text{min}}\_l) \right)^{\frac{1}{2}} - q^2 \left( \sum\_{l=1}^L \mathsf{Tr}(\mathsf{W}\_l) + P\_0 \right) + \lambda\_1 \sum\_{l=1}^L \mathsf{P}(\eta\_{k,l}), \tag{51}$$

$$s.t. \ C\_2: \Gamma\_{k,l}(\mathbf{W}) - \mathbf{I}\_{k,l}(\mathbf{W}\_\prime \mathcal{L}\_{k,l}) \le \mathbf{0}\_\prime \,\forall k, l\_\prime \tag{52}$$

$$\text{C}\_{3\text{\textquotedblleft C}4\text{\textquotedblright}C\_5}\text{C}\_6\text{\textquotedblleft C}\_7\text{\textquotedblright}\text{\textquotedblleft}Q\_{4\text{\textquotedblleft}Q\_{5\text{\textquotedblright}}}\tag{53}$$

$$\mathbf{C}\_{8}: 4\zeta\_{k,l} + \left(\eta\_{k,l} - \text{SINR}\_{l}^{\text{min}}\right)^{2} \ge \left(\eta\_{k,l} + \text{SINR}\_{l}^{\text{min}}\right)^{2} \forall k, l,\tag{54}$$

#### *4.3. CCCP Algorithm*

From (51), (52) and (54), it can be seen that the problem *Q*<sup>5</sup> is a DC programming problem. To handle this challenge, the CCCP framework algorithm is a common method to solve the DC programming problem [35], which is an iterative framework including two operations: convexification and optimization. In the convexification step, by adopting the first-order Taylor expansion, the convex part of the objective function and the concave part of the constraint function can be linearized; then, the DC programming problem is transformed into a convex problem. It should be noted that the convex problem obtained from the convexification step provides a global lower bound for the original problem, and the optimization step is mainly to maximize the lower bound. Meanwhile, the performance of the CCCP algorithm is closely related to the initial point of the variables, but the equality constraint *C*<sup>4</sup> of the problem *Q*<sup>5</sup> limits the selection of the initial point. To find a feasible initial point, we substitute the constraint *C*<sup>4</sup> into the objective function and set the penalty factor *λ*<sup>2</sup> > 0. Then, the problem *Q*<sup>5</sup> can be equivalently converted to:

$$\begin{split} Q\_6: \max\_{\eta, \mathbf{W} \text{SINR}^{\text{min}}\_l, \mathbf{J}} & \text{EE} = 2q \left( B \sum\_{l=1}^L \mathbf{Y}\_l (\text{SINR}^{\text{min}}\_l) \right)^{\frac{1}{2}} - q^2 \left( \sum\_{l=1}^L \text{Tr} (\mathbf{W}\_l) + P\_0 \right) \\ & + \lambda\_1 \sum\_{l=1}^L \mathbf{P} (\eta\_{k,l}) - \sum\_{l}^L \lambda\_2 \left( \sum\_{k}^{|\mathcal{U}\_l|} \eta\_{k,l} - \mathcal{U}\_\mathbf{s} \right)^2 \end{split} \tag{55}$$

$$\text{s.t.}\,\mathsf{C}\_{2}\,\mathsf{C}\_{3}\,\mathsf{C}\_{5}\,\mathsf{C}\_{6}\,\mathsf{C}\_{7}\,\mathsf{C}\_{8}\,\text{in}\,\mathrm{Q}\_{5},\tag{56}$$

*Convexification*: Let *ηk*,*<sup>l</sup>* , SINRmin *l* ,*W*, *ζ* (*t*−1) be the estimated value of variables *ηk*,*<sup>l</sup>* , SINRmin *l* ,*W*, *ζ* in iteration *t* − 1 of the problem *Q*6. In iteration *t*, for the convex part *λ*1 *L* ∑ *l*=1 P(*ηk*,*l*), we adopt the first-order Taylor expansion to replace it, which is reflected in line three of Algorithm 1. The first-order Taylor expansion of P(*ηk*,*l*) can be expressed as:

$$\mathbf{P}(\eta\_{k,l})^{te} = \mathbf{P}\left(\eta\_{k,l}^{(t-1)}\right) + \left(\eta\_{k,l} - \eta\_{k,l}^{(t-1)}\right)\nabla\mathbf{P}\left(\eta\_{k,l}^{(t-1)}\right),\tag{57}$$

$$\nabla \mathcal{P} \left( \eta\_{k,l}^{(t-1)} \right) = \log \left( \eta\_{k,l}^{(t-1)} \right) - \log \left( 1 - \eta\_{k,l}^{(t-1)} \right),\tag{58}$$

Similarly, in the constraint *C*<sup>2</sup> : Γ*k*,*l*(*W*) − I*k*,*l*(*W*, *ζk*,*l*) ≤ 0, we replace the concave function I*k*,*l*(*W*, *ζk*,*l*) with its first-order Taylor expansion, which is reflected in line three of Algorithm 1. The first-order Taylor expansion of I*k*,*l*(*W*, *ζk*,*l*) can be expressed as:

$$\mathbf{I}\_{k,l}(\mathbf{W}, \boldsymbol{\zeta}\_{k,l})^{t\varepsilon} = \mathbf{I}\_{k,l}\left(\mathbf{W}^{(t-1)}, \boldsymbol{\zeta}\_{k,l}^{(t-1)}\right) + \nabla^{T}\mathbf{I}\_{k,l}\left(\mathbf{W}^{(t-1)}, \boldsymbol{\zeta}\_{k,l}^{(t-1)}\right) \begin{bmatrix} \left\{\mathbf{W}\_{l} - \mathbf{W}\_{l}^{(t-1)}\right\}\_{l=1}^{L} \\ \boldsymbol{\zeta}\_{k,l} - \boldsymbol{\zeta}\_{k,l}^{(t-1)} \end{bmatrix},\tag{59}$$

$$\nabla \mathbf{I}\_{k,l}\left(\mathbf{W}^{(t-1)}, \boldsymbol{\zeta}\_{k,l}^{(t-1)}\right) = \left[ \left\{ \frac{\left(\overline{\mathbf{H}}\_{k,l}\right)^{T}}{1 + \boldsymbol{\zeta}\_{k,l}^{(t-1)}} \right\}\_{l=1}^{L} - \frac{\sigma^{2} + \sum\_{l=1}^{L} \text{Tr}\left(\overline{\mathbf{H}}\_{k,l} \mathbf{W}\_{l}^{(t-1)}\right)}{\left(1 + \boldsymbol{\zeta}\_{k,l}^{(t-1)}\right)^{2}} \right]^{T},\tag{60}$$

In the constraint *C*<sup>8</sup> : 4*ζk*,*<sup>l</sup>* + *<sup>η</sup>k*,*<sup>l</sup>* <sup>−</sup> SINRmin *l* 2 ≥ *ηk*,*<sup>l</sup>* + SINRmin *l* 2 ∀*k*, *l*, we replace the concave function *<sup>η</sup>k*,*<sup>l</sup>* <sup>−</sup> SINRmin *l* 2 with its first-order Taylor expansion, which is reflected in line 3 of Algorithm 1. The first-order Taylor expansion of *<sup>η</sup>k*,*<sup>l</sup>* <sup>−</sup> SINRmin *l* 2 can be expressed as:

$$\begin{aligned} \left(\eta\_{k,l} - \text{SINR}\_{l}^{\text{min}}\right)^{2,t\varepsilon} &= \left(\eta\_{k,l}^{(t-1)} - \text{SINR}\_{l}^{\text{min},(t-1)}\right)^{2} + \begin{bmatrix} 2\left(\eta\_{k,l}^{(t-1)} - \text{SINR}\_{l}^{\text{min},(t-1)}\right) \\ -2\left(\eta\_{k,l}^{(t-1)} - \text{SINR}\_{l}^{\text{min},(t-1)}\right) \end{bmatrix}^{T} \begin{bmatrix} \eta\_{k,l} - \eta\_{k,l}^{(t-1)} \\ \text{SINR}\_{l}^{\text{min}} - \text{SINR}\_{l}^{\text{min},(t-1)} \end{bmatrix}, \end{aligned} \tag{61}$$

*Optimization*: The optimization step is reflected in line nine of Algorithm 1. According to (57), (59) and (61), the problem *Q*<sup>6</sup> can be equivalently converted to

$$\begin{split} Q\_7: \max\_{\eta, \mathbf{W}, \text{SINR}^{\text{min}}\_l, \mathcal{J}} \text{EE} &= \quad 2q \left( B \sum\_{l=1}^L \mathbf{Y}\_l (\text{SINR}^{\text{min}}\_l) \right)^{\frac{1}{2}} - q^2 \left( \sum\_{l=1}^L \text{Tr} (\mathbf{W}\_l) + \mathbf{P}\_0 \right) \\ &+ \lambda\_1 \sum\_{l=1}^L \mathbf{P} (\eta\_{k,l})^{t\varepsilon} - \sum\_{l}^L \lambda\_2 \left( \sum\_{k}^{|\mathcal{U}\_l|} \eta\_{k,l} - \mathcal{U}\_s \right)^2 \end{split} \tag{62}$$

*s*.*t*. *C*3, *C*<sup>5</sup> ,*C*<sup>6</sup> , *C*<sup>7</sup> *in Q*6, (63)

$$\mathbf{C}\_{2}: \Gamma\_{k,l}(\mathbf{W}) - \mathbf{I}\_{k,l}(\mathbf{W}, \mathbb{Q}\_{k,l})^{t\varepsilon} \le \mathbf{0}, \forall k, l,\tag{64}$$

$$\mathbf{C}\_{8}: 4\zeta\_{k,l} + \left(\eta\_{k,l} - \text{SINR}\_{l}^{\text{min}}\right)^{2,te} \ge \left(\eta\_{k,l} + \text{SINR}\_{l}^{\text{min}}\right)^{2} \forall k, l,\tag{65}$$

For the problem *<sup>Q</sup>*7, the variables *ηk*,*<sup>l</sup>* , SINRmin *l* ,*W*, *ζ* (*t*+1) can be updated by the iterative optimization.

*Feasible Initial Point*: It should be noted that the CCCP algorithm needs a feasible initial point to ensure that the algorithm converges to a stationary point, as the selection of the initial point can affect the performance of the CCCP algorithm. To find a better initial point, we adopt the following method, which is reflected in line one of Algorithm 1.


$$P\_{FSS} : \left\{ \mathbf{W}^{(0)} \right\} : \min\_{\mathbf{W}} \sum\_{l=1}^{L} \text{Tr}(\mathbf{W}\_l) \,\tag{66}$$

$$\left\| s.t. \mathcal{C}\_1 : \left\| \sigma^2 \dots .Tr(\overline{\mathbf{H}}\_{k,l} \mathbf{W}\_j)\_{j \neq l} \dots \right\| \le \frac{Tr(\overline{\mathbf{H}}\_{k,l} \mathbf{W}\_l)}{\eta\_{k,l} \text{SINR}\_0} , \forall k, l \right\} \tag{67}$$

;

$$\mathbf{C}\_{2}: \sum\_{l=1}^{L} \operatorname{Tr}(\mathbf{W}\_{l}) \le P\_{T\prime} \tag{68}$$


and update  $\eta\_{k,l}^{(0)}$  according to  $\eta\_{k,l}^{(0)} = \min\left\{1, \frac{\text{SINR}\_{k,l}^{(0)}}{\text{SINR}\_{0}}\right\}$ 

• Based on *η* (0) *k*,*l* and *W*(0) , calculate *ζ* (0) and <sup>n</sup> SINRmin, *l* (0) o*<sup>L</sup> l*=1 .

#### *4.4. Penalty Iteration Algorithm*

It should be noted that the SDP algorithm brings the nonconvex and nonsmooth constraint, i.e., *C*<sup>7</sup> : *rank*(*W<sup>l</sup>* ) = 1. To solve the rank-one constraint, many existing research directly relaxes the rank-one constraint in the optimization step [36], and then judges whether the optimization solution {*Wl*} *L <sup>l</sup>*=<sup>1</sup> meets the rank-one constraint. If so, the eigenvalue decomposition (EVD) algorithm is directly adopted to obtain the hybrid beamforming vectors {*f l* } *L l*=1 according to *W<sup>l</sup>* = *f l f H l* , and if the optimization solution {*Wl*} *L l*=1 does not meet the rank-one constraint, the Gaussian randomization algorithm (GRA) is usually adopted. The basic idea of the GRA is as follows: Firstly, a set of candidate Gaussian vectors n *wg*,*<sup>l</sup>* o*<sup>G</sup> g*=1 *<sup>L</sup> l*=1 are generated based on the optimization solution

{*Wl*} *L l*=1 , where *G* represents the number of the Gaussian randomization. Secondly, from the generated G-group candidate Gaussian vector pool, combined with the power redistribution among the multicast groups, a group of Gaussian vectors is selected as the optimal hybrid beamforming matrix to maximize the objective function in the problem *Q*7. It should be noted that in the case of the high-dimensional matrix, GRA has high complexity and large performance loss, resulting in poor availability. To this end, we adopt a feasible algorithm with the better performance: the penalty iteration algorithm.

According to the properties of the matrix, *rank*(*W<sup>l</sup>* ) = 1 is equivalent to *Tr*(*W<sup>l</sup>* ) − *λ*max(*Wl*) = 0. Therefore, the nonsmooth method is adopted to transform the constraint *rank*(*W<sup>l</sup>* ) = 1 in the problem *Q*<sup>7</sup> into the following form:

$$C\_7: \operatorname{Tr}(\mathbf{W}\_l) - \lambda\_{\max}(\mathbf{W}\_l) \le 0,\tag{69}$$

where *λ*max(*Wl*) is the function of solving the maximum eigenvalue. It should be noted that for any positive semidefinite matrix *W<sup>l</sup>*0, the inequality *Tr*(*W<sup>l</sup>* ) − *λ*max(*Wl*) ≥ 0 is always true, which means that the transformed constraint *C*<sup>7</sup> and *Tr*(*W<sup>l</sup>* ) − *λ*max(*Wl*) = 0 are equivalent. Then, we can obtain that the matrix *W<sup>l</sup>* has only one non-zero eigenvalue and can be given by

$$\mathbf{W}\_{l} = \lambda\_{\text{max}}(\mathbf{W}\_{l}) \mathbf{w}\_{l,\text{max}} \mathbf{w}\_{l,\text{max}}^{H} \tag{70}$$

where *wl*,max is the corresponding unit eigenvector. Therefore, the problem *Q*<sup>7</sup> can be converted into

$$\begin{split} \text{Qs}: \max\_{\eta, \mathbf{W}, \text{SINR}^{\text{min}}\_{l}, \mathcal{L}} \text{EE} &= \quad 2q \left( B \sum\_{l=1}^{L} \mathbf{Y}\_{l} (\text{SINR}^{\text{min}}\_{l}) \right)^{\frac{1}{2}} - q^{2} \left( \sum\_{l=1}^{L} \text{Tr} (\mathbf{W}\_{l}) + \mathbf{P}\_{0} \right) \\ &+ \lambda\_{1} \sum\_{l=1}^{L} \text{P} (\eta\_{k,l})^{\text{te}} - \sum\_{l}^{L} \lambda\_{2} \left( \sum\_{k}^{|\mathcal{U}\_{l}|} \eta\_{k,l} - \mathcal{U}\_{\mathbf{s}} \right)^{2} \end{split} \tag{71}$$

$$\text{is.t.}\,\mathsf{C}\_{2}\,\mathsf{C}\_{3}\,\mathsf{C}\_{5}\,\mathsf{C}\_{6}\,\mathsf{C}\_{6}\,\mathsf{C}\_{8}\,\text{in}\,\mathsf{Q}\_{7},\tag{72}$$

$$\mathcal{C}\_{7}: Tr(\mathbf{W}\_{l}) - \lambda\_{\text{max}}(\mathbf{W}\_{l}) \le 0,\tag{73}$$

In the iterative calculation of the problem *Q*8, based on the obtained *W<sup>l</sup>* , if the value of *Tr*(*W<sup>l</sup>* ) − *λ*max(*Wl*) is small enough, the matrix *W<sup>l</sup>* can be considered to meet the rank-one constraint, which is reflected in line seven of Algorithm 1. Therefore, to make the value of *Tr*(*W<sup>l</sup>* ) − *λ*max(*Wl*) as small as possible, we adopt the penalty iteration algorithm and substitute the constraint *C*<sup>7</sup> into the objective function in the problem *Q*8. Therefore, the problem *Q*<sup>8</sup> can be converted into

$$Q\_{\mathcal{Y}}: \max\_{\boldsymbol{\eta}, \mathbf{W}, \text{SINR}^{\text{min}}\_{l}, \boldsymbol{\zeta}} \text{EE} = F\left(\boldsymbol{\eta}, \mathbf{W}, \text{SINR}^{\text{min}}\_{l}, \boldsymbol{\zeta}\right) - \lambda\_3 \sum\_{l}^{L} (\text{Tr}(\mathbf{W}\_l) - \lambda\_{\text{max}}(\mathbf{W}\_l))\_\ast \tag{74}$$

$$\text{s.t.}\ \mathsf{C}\_{2}\ \mathsf{C}\_{3}\ \mathsf{C}\_{5}\ \mathsf{C}\_{6}\ \mathsf{C}\_{6}\ \mathsf{C}\_{8}\ in\ \mathsf{Q}\_{8}\tag{75}$$

$$\begin{array}{rcl} \text{where} & \begin{array}{rcl} F\left(\boldsymbol{\eta}, \boldsymbol{\mathsf{W}}, \boldsymbol{\mathsf{S}} \text{Sign} \boldsymbol{\mathsf{R}}\_{l}^{\min}, \boldsymbol{\zeta} \right) & = & \mathbf{2}q \left( \boldsymbol{B} \sum\_{l=1}^{L} \mathbf{Y}\_{l} (\boldsymbol{\mathsf{S}} \text{Sign} \boldsymbol{\mathsf{R}}\_{l}^{\min}) \right)^{\frac{1}{2}} - q^{2} \left( \sum\_{l=1}^{L} \text{Tr} (\boldsymbol{\mathsf{W}}\_{l}) + \boldsymbol{\mathsf{P}}\_{0} \right) \\ + \boldsymbol{\lambda} \cdot \boldsymbol{\lambda} \cdot \sum^{L} \mathbf{P} (\boldsymbol{\mathsf{w}}, \boldsymbol{\lambda})^{\text{fl}} - \sum^{L} \boldsymbol{\lambda} \cdot \begin{pmatrix} |\boldsymbol{\mathcal{U}}| & & \boldsymbol{\mathsf{W}} \end{pmatrix}^{2} \end{array} \text{ and } \boldsymbol{\lambda} > 0 \text{ is the penalty factor, which is constant.}$$

+ *λ*<sup>1</sup> ∑ *l*=1 P(*ηk*,*l*) *te* <sup>−</sup> ∑ *l λ*2 ∑ *k ηk*,*<sup>l</sup>* − *U<sup>s</sup>* and *λ*<sup>3</sup> > 0 is the penalty factor, which is gen-

erally larger enough to ensure that a smaller value of *Tr*(*W<sup>l</sup>* ) − *λ*max(*Wl*) can be obtained. According to (74), the iterative calculation of the problem *Q*<sup>9</sup> can maximize function

*F η*,*W*, SINRmin *l* , *ζ* and minimize function *Tr*(*W<sup>l</sup>* ) − *λ*max(*Wl*). It should be noted that *Tr*(*W<sup>l</sup>* ) is an affine function and *λ*max(*Wl*) is nonsmooth, which can result in the nonsmoothness of the objective function in the problem *Q*9. To handle the challenge, we replace *λ*max(*Wl*) with its first-order Taylor expansion. The subgradient of *λ*max(*Wl*) is *∂λ*max(*W<sup>l</sup>* ) *∂W<sup>l</sup>* = *wl*,max*w<sup>H</sup> <sup>l</sup>*,max and its first-order Taylor expansion can be expressed as follows, which is reflected in line eight of Algorithm 1.

$$
\lambda\_{\max} \left( \mathbf{W}\_l^{(t)} \right) \ge \lambda\_{\max} \left( \mathbf{W}\_l^{(t-1)} \right) + \left< \mathbf{w}\_{l,\max} \mathbf{w}\_{l,\max}^H \mathbf{W}\_l^{(t)} - \mathbf{W}\_l^{(t-1)} \right>,\tag{76}
$$

$$\text{where}\left\langle \mathfrak{w}\_{l,\text{max}} \mathfrak{w}\_{l,\text{max}}^H \mathbf{W}\_l^{(t)} - \mathbf{W}\_l^{(t-1)} \right\rangle = \text{Tr}\left( \left( \mathfrak{w}\_{l,\text{max}} \mathfrak{w}\_{l,\text{max}}^H \right)^H \left( \mathbf{W}\_l^{(t)} - \mathbf{W}\_l^{(t-1)} \right) \right).$$

We substitute (76) into the objective function in the problem *Q*<sup>9</sup> to replace *λ*max(*Wl*), and the problem *Q*<sup>9</sup> can be expressed as:

$$\mathbf{Q}\_{10}: \max\_{\boldsymbol{\eta}, \mathbf{W}, \text{SNR}\_{l}^{\text{min}}, \boldsymbol{\zeta}} \operatorname{EE} = \mathbf{F} \left( \boldsymbol{\eta}, \mathbf{W}, \text{SNR}\_{l}^{\text{min}}, \boldsymbol{\zeta} \right) - \lambda\_3 \sum\_{l}^{L} \left( \operatorname{Tr} (\mathbf{W}\_{l}) - \lambda\_{\max} \left( \mathbf{W}\_{l}^{(t-1)} \right) + \left< \mathbf{w}\_{l, \max} \mathbf{w}\_{l, \max}^{H}, \mathbf{W}\_{l}^{(t)} - \mathbf{W}\_{l}^{(t-1)} \right> \right), \tag{77}$$

*s*.*t*. *C*2, *C*3, *C*<sup>5</sup> ,*C*<sup>6</sup> , *C*<sup>8</sup> *in Q*9, (78)

In conclusion, the robust joint user scheduling and hybrid beamforming design algorithm for the massive MIMO LEO satellite multigroup multicast communication system is shown in Algorithm 1.

**Algorithm 1:** Joint user scheduling and hybrid beamforming design algorithm.

.

**Input:** CCCP algorithm iteration index *k*, thresholds *ε*<sup>1</sup> , penalty iteration algorithm iteration index *m*, thresholds *ε*2, penalty factor *λ*<sup>1</sup> , *λ*2, *λ*3.

$$\text{1. Initial: } \left(\mathfrak{h}\_{\prime}\mathbf{W}\_{\prime}\mathbf{S}\mathbf{I}\mathbf{N}\mathbf{R}\_{I}^{\min}, \mathfrak{f}\right)^{(k=0)}, q^{(k=0)}$$


8. Calculate the maximum eigenvalue *λ*max(*W<sup>l</sup>* ) of *W* (*m*) *l* and the corresponding eigenvector *w* (*m*) *l*,max .

.

9. Using CVX toolbox, calculate the variables *η*,*W*, SINRmin *l* , *ζ* (*m*) *opt* at the *m*th iteration according to (77).

10. If <sup>n</sup> *W* (*m*+1) *l* o*<sup>L</sup> l* ≈ n *W* (*m*) *l* o*<sup>L</sup> l* , then 11. Update *λ*<sup>3</sup> = 2*λ*3. 12. else 13. Update *m* = *m* + 1. 14. end 15. end 16. Update *η*,*W*, SINRmin *l* , *ζ* (*k*+1) = *η*,*W*, SINRmin *l* , *ζ* (*m*) , *k* = *k* + 1, *λ*<sup>1</sup> = *λ*<sup>1</sup> + 1, *λ*<sup>2</sup> = *λ*<sup>2</sup> + 1. 17. end **Output**: *η*, *W*, *SINR*min *l* , *ζ opt*.

#### **5. Convergence and Complexity Analysis**

#### *5.1. Convergence*

The effectiveness of the algorithm depends on its convergence. For the convergence of the CCCP algorithm, the convergence has been proven by [37]. To prove the convergence of the penalty iteration algorithm, let the variable solution and objective function value of the optimization problem *<sup>Q</sup>*<sup>10</sup> be *η*,*W*, SINRmin *l* , *ζ* (*k*+1) and *<sup>F</sup>EE η*,*W*, SINRmin *l* , *ζ* (*k*+1) at the *k*th iteration. Therefore, the convergence can be proved as follows:

$$\begin{split} &F\_{\text{EE}}\left(\left(\boldsymbol{\eta},\mathbf{W},\text{SINR}\_{l}^{\text{min}},\boldsymbol{\zeta}\right)^{(k+1)}\right) = F\left(\left(\boldsymbol{\eta},\mathbf{W},\text{SINR}\_{l}^{\text{min}},\boldsymbol{\zeta}\right)^{(k+1)}\right) - \lambda\_{3}\frac{L}{l}\Big{(}\text{Tr}\left(\mathbf{W}\_{l}^{(k+1)}\right) - \lambda\_{\text{max}}\left(\mathbf{W}\_{l}^{(k+1)}\right)\Big{)} \\ &\geq F\left(\left(\boldsymbol{\eta},\mathbf{W},\text{SINR}\_{l}^{\text{min}},\boldsymbol{\zeta}\right)^{(k+1)}\right) - \lambda\_{3}\frac{L}{l}\Big{(}\text{Tr}\left(\mathbf{W}\_{l}^{(k+1)}\right) - \lambda\_{\text{max}}\left(\mathbf{W}\_{l}^{(k)}\right) - \left(\mathbf{w}\_{l,\text{max}}\mathbf{w}\_{l,\text{max}}^{H},\mathbf{W}\_{l}^{(k+1)} - \mathbf{W}\_{l}^{(k)}\right)\Big{)} \\ &\overset{\text{by}(7)}{\geq} F\left(\left(\boldsymbol{\eta},\mathbf{W},\text{SINR}\_{l}^{\text{min}},\boldsymbol{\zeta}\right)^{(k)}\right) - \lambda\_{3}\frac{L}{l}\Big{(}\text{Tr}\left(\mathbf{W}\_{l}^{(k)}\right) - \lambda\_{\text{max}}\left(\mathbf{W}\_{l}^{(k)}\right)\Big{)} \\ &= F\_{\text{EE}}\left(\left(\boldsymbol{\eta},\mathbf{W},\text{SINR}\_{l}^{\text{min}},\boldsymbol{\zeta}\right)^{(k)}\right) \end{split} \tag{79}$$

The convergence can be proved according to (79). Therefore, after initializing the values of *η*,*W*, SINRmin *l* , *ζ* (*k*=0) , *λ* (*k*=0) 1 , *λ* (*k*=0) 2 and *λ* (*k*=0) 3 , the proposed algorithm can iteratively converge to an optimal solution by setting a reasonable convergence threshold.

#### *5.2. Complexity*

The complexity of the algorithm directly affects its performance. In the algorithms adopted, the complexity of the hierarchical clustering algorithm can be calculated according to the connection algorithm, similarity measurement criteria and hierarchical grouping process, and the algorithm complexity can be expressed as *O LK*2*N* . The complexity of the joint user scheduling and hybrid beamforming design algorithm is closely related to the number of multicast groups and scheduling users. In addition, the number of optimization variables and constraints in the CCCP algorithm and the penalty iteration algorithm can also affect the complexity [38]. In the problem *Q*10, the number of optimization variables is 2 *L* ∑ *l*=1 |*Ul* | + 2*L*, the number of convex constraints is *L* ∑ *l*=1 |*Ul* | and the number of linear constraints is *L* ∑ *l*=1 |*Ul* | + 3*L*. Let the number of iterations in the penalty iteration algorithm and the CCCP algorithm be *I<sup>p</sup>* and *Ic*, respectively. In conclusion, the overall complexity of the proposed algorithm is *O LK*2*N* + *I<sup>p</sup> I<sup>c</sup>* 2 *L* ∑ |*Ul* | + 2*L* <sup>2</sup> *L* ∑ |*Ul* | + 3*L* .

*l*=1 *l*=1 According to the algorithm complexity, the proposed joint user scheduling and hybrid beamforming design algorithm has a strong timeliness in small dimensional communication systems. However, for large dimensional communication systems, such as the satellite communication system, the number of active users is usually large. According to the complexity analysis, with the increase in the total number of active users, the convergence speed of the algorithm would gradually slow down, and the complexity would gradually increase. Considering the characteristics of the LEO satellite communication system, the delay caused by the high complexity is unacceptable, which would affect the overall performance of the communication system. To handle this problem, considering the balance of the algorithm performance and complexity, before the joint user scheduling and hybrid beamforming design, we can appropriately reduce the system dimension by adding the user preselection step in the algorithm process. The user preselection step can reduce the total number of active users in each transmission, and then we carry out the joint user scheduling and hybrid beamforming design for preselected users.

#### **6. User Preselection Algorithm**

In the user preselection step, *Ul*,*<sup>p</sup>* users in the *l*th multicast group are preselected as the user representatives, where *U<sup>s</sup>* < *Ul*,*<sup>p</sup>* ≤ |*U<sup>l</sup>* |. The user preselection process is shown in Figure 6. The symbols (circles, squares and triangles) represent the users in different the multicast group. The red circles represent preselected users and scheduling users.

*5.2. Complexity*

variables is

linear constraints is

1 2 2 *L l l*

**6. User Preselection Algorithm**

as the user representatives, where

In the user preselection step,

rithm and the CCCP algorithm be

complexity of the proposed algorithm is

= <sup>+</sup>

*U L*

1

= <sup>+</sup>

*L l l*

**Figure 6.** Design process of joint user scheduling and beamforming with the user preselection. **Figure 6.** Design process of joint user scheduling and beamforming with the user preselection.

The complexity of the algorithm directly affects its performance. In the algorithms adopted, the complexity of the hierarchical clustering algorithm can be calculated according to the connection algorithm, similarity measurement criteria and hierarchical group-

of the joint user scheduling and hybrid beamforming design algorithm is closely related to the number of multicast groups and scheduling users. In addition, the number of optimization variables and constraints in the CCCP algorithm and the penalty iteration algo-

, the number of convex constraints is

and *c I*

2

According to the algorithm complexity, the proposed joint user scheduling and hybrid beamforming design algorithm has a strong timeliness in small dimensional communication systems. However, for large dimensional communication systems, such as the satellite communication system, the number of active users is usually large. According to the complexity analysis, with the increase in the total number of active users, the convergence speed of the algorithm would gradually slow down, and the complexity would gradually increase. Considering the characteristics of the LEO satellite communication system, the delay caused by the high complexity is unacceptable, which would affect the overall performance of the communication system. To handle this problem, considering the balance of the algorithm performance and complexity, before the joint user scheduling and hybrid beamforming design, we can appropriately reduce the system dimension by adding the user preselection step in the algorithm process. The user preselection step can reduce the total number of active users in each transmission, and then we carry out the

*p I*

joint user scheduling and hybrid beamforming design for preselected users.

*Ul p*,

*U U U s l p l* ,

users in the

Figure 6. The symbols (circles, squares and triangles) represent the users in different the multicast group. The red circles represent preselected users and scheduling users.

*l*th

( ) <sup>2</sup> *O LK N*

1 *L l l U* = 

, respectively. In conclusion, the overall

.

multicast group are preselected

. The user preselection process is shown in

1 1 2 2 2 3 *L L p c <sup>l</sup> <sup>l</sup> l l O LK N I I U L U L* = = <sup>+</sup> <sup>+</sup> <sup>+</sup> 

*Q*<sup>10</sup>

. Let the number of iterations in the penalty iteration algo-

. The complexity

and the number of

, the number of optimization

ing process, and the algorithm complexity can be expressed as

rithm can also affect the complexity [38]. In the problem

3

*U L*

The selection of preselected users can affect the performance of the joint user scheduling and hybrid beamforming design algorithm, which depends on the preselected algorithm. In the beamforming design of the multigroup multicast communication system, the beamforming vector is oriented to multiple users in the multicast group. Therefore, in the process of user preselection, to maximize the receive gain of each user, i.e., *h H k*,*l f l* , the beamforming vector *f l* of the *l*th multicast group should be collinear with the users' channel vectors in the multicast group as far as possible. Therefore, the channel vectors of the preselected users in the same multicast group should also be strongly linearly correlated. Meanwhile, the interference among multicast groups should also be taken into account in the user preselection stage. To reduce the interference among multicast groups, the channel vectors of preselected users among different multicast groups should be orthogonal. Similarly, the beamforming vector of the multicast group should be orthogonal to the users' channel vectors in other multicast groups.

In conclusion, we adopt a low complexity user preselection algorithm, which can preselect orthogonal users among the different multicast groups and linearly correlated users in the same multicast group. The proposed algorithm is divided into two steps, as follows:


The specific preselection process of the two steps is as follows:

**Algorithm 2:** User preselection algorithm.

**Step 1: Orthogonal user preselection algorithm among the different multicast groups. Input**: CSI.

1. Let Id(1) = *Index* max *hk*,*<sup>l</sup>* , ∀*k*, *l* , select the user with the largest channel gain, Id(1) is the index of the user.

2. while *l* ≤ *L*, *l* 6= Id(1)

3. For all users in the *l*th multicast group, calculate *Zk*,*<sup>l</sup>* = *hk*,*<sup>l</sup> I<sup>N</sup>* − ∑ Id(*l*) *j*=Id(1) *h H* (*j*) *h*(*j*) k*h*(*j*)k 2 in turn.

2 4. Id(*l*) = *Index* max *Zk*,*<sup>l</sup>* , ∀*k* ∈ *l* , the user with index Id(*l*) is the preselected orthogonal user of the *l*th multicast group.

5. end

**Output**: Orthogonal users among the different multicast groups.

**Step 2: User preselection algorithm in each multicast group.**

**Input**: Orthogonal users among the different multicast groups, CSI.

1. For *l* = 1 : *L*

2. For other users in the *l*th multicast group except the orthogonal user preselected in step 1, calculate the linear correlation value between each user and the preselected orthogonal user of the

multicast group in turn, i.e., *Ck*,*<sup>l</sup>* = *h H k*,*l <sup>h</sup>*Id(*l*) *h H* Id(*l*) *h H* Id(*l*) 2 2 .

3. Based on the *<sup>C</sup>k*,*<sup>l</sup>* of users in each multicast group, select top *Ul*,*<sup>p</sup>* − 1 largest users, plus the orthogonal users in step 1 as the preselected users of each multicast group.

4. end

5. end

**Output**: Preselected users for each multicast group.

After the user preselection, the joint user scheduling and hybrid beamforming design is for the preselected users. Therefore, the dimension of the LEO satellite communication system will be reduced, and the algorithm complexity will be reduced. Although the algorithm performance has a slight loss, compared with the decoupling design of user scheduling and beamforming, the performance is greatly improved. In conclusion, the joint user scheduling and hybrid beamforming design with the user preselection step is a better choice after balancing performance and complexity.

#### **7. Solution of The Digital Beamforming Matrix and The Analog Beamforming Matrix**

In this section, we aim to investigate the design of digital beamforming matrix *FBB* and analog beamforming matrix *FRF* in a hybrid beamformer. After obtaining *W*, we need to further solve *FBB* and *FRF*. The solution method of *FBB* and *FRF* can be divided into two steps:


#### *7.1. Solution of The Hybrid Beamforming Matrix*

Before calculating *FBB* and *FRF*, it is necessary to obtain the hybrid beamforming matrix *F*. For the solution of *F*, we can adopt the EVD algorithm based on the previously obtained optimization variable *W*. According to the relationship between *W* and *F*, i.e., *W<sup>l</sup>* , *f l f H l* , the solution of *F* can be modeled as follows:

$$\min\_{f\_l} \left\| \mathbf{W}\_l - f\_l f\_l^H \right\|\_{F'}^2 \tag{80}$$

where the hybrid beamforming vector *f l* can be given by

$$f\_l = \sqrt{v\_l} \mathfrak{u}\_{l\prime} \,\forall l\prime \,\tag{81}$$

where *v<sup>l</sup>* is the maximum eigenvalue of the matrix *W<sup>l</sup>* and *u<sup>l</sup>* is the maximum eigenvector of the matrix *W<sup>l</sup>* .

## *7.2. MM-AltOpt Algorithm: Solution of FBB and FRF*

According to *F* = *FRFFBB* and (*FRF*)*i*,*<sup>j</sup>* <sup>=</sup> 1, the solution of *<sup>F</sup>BB* and *<sup>F</sup>BB* can be modeled as a joint optimization problem with the power and constant modulus constraints, as follows:

$$P\_1: \min\_{\mathbf{F}\_{BB}, \mathbf{F}\_{RF}} \left\| \mathbf{F}^{opt} - \mathbf{F}\_{RF} \mathbf{F}\_{BB} \right\|\_{\mathbf{F}'}^2 \tag{82}$$

$$\text{s.t.}\,\mathsf{C}\_{1}:\mathsf{F}\_{RF}\in\mathbb{F}\_{\prime}\tag{83}$$

$$\mathbf{C}\_{2}: \|\mathbf{F}\_{\text{RF}}\mathbf{F}\_{\text{BB}}\|\_{F}^{2} = P\_{T\prime} \tag{84}$$

where F = *FRF* ∈ *C N*×*L FRF*(*i*,*j*) 2 = 1, 1 ≤ *i* ≤ *N*, 1 ≤ *j* ≤ *L* represents the unit modulus constraint, which is determined by the phase shifter in UPA, and the constraint *C*<sup>2</sup> represents the power constraint.

It is worth noting that the problem *P*<sup>1</sup> is a matrix decomposition problem with the constant modulus constraint and the equality constraint. The objective function is a nonconvex function of variables *FBB* and *FRF*, and the constraints *C*<sup>1</sup> and *C*<sup>2</sup> are also nonconvex. Meanwhile, it can be seen that when one of the two variables is given, the objective function is the convex function of the other variable. To solve the problem *P*1, we invoke the alternating optimization algorithm. The alternating optimization algorithm can decompose the multivariable joint optimization problem into multiple subproblems according to the partial convexity of the problem *P*1, and one of the variables can be iteratively solved by fixing the residual variables.

It should be noted that the nonconvexity of constraints is still a challenge. To this end, we first relax the constraint *C*2, and then use the scale factor to adjust the digital beamforming matrix *FBB* to meet the power constraint. Then, for the solution of the analog beamforming matrix *FRF* with the unit modulus constraint, the MM algorithm is adopted [33].

#### *7.3. Solution of The Analog Beamforming Matrix Based on The MM Algorithm*

According to the solution process of the alternating optimization algorithm, we first solve the analog beamforming matrix *FRF* based on the digital beamforming matrix *FBB*. Thus, the problem *P*<sup>1</sup> can be expressed as:

$$P\_2: \min\_{F\_{RF}} \left\| F - F\_{RF} F\_{BB}^{(n)} \right\|\_{F}^{2} \tag{85}$$

$$\text{s.t.}\,\mathsf{C}\_{1}:\mathsf{F}\_{\mathsf{RF}}\in\mathbb{F}\_{\mathsf{A}}\tag{86}$$

where *F* (*n*) *BB* represents the estimated value of the digital beamforming matrix *FBB* at the *n*th iteration. Due to the unit module constraint of elements in *FRF*, the problem *P*<sup>2</sup> is a nonconvex optimization problem.

According to the MM framework theory, the key step is constructing a surrogate function of the objective function in the optimization problem [39]. To construct the surrogate function, we decompose the matrix *F* by rows. According to the equivalence of the F-norm and *L*2-norm of the vector, the problem *P*<sup>2</sup> can be rewritten as:

$$P\_3: \min\_{\mathbf{F}\_{\rm RF}} \sum\_{i=1}^{N} \mathbf{F}\_i^H \mathbf{F}\_i - 2\Re\left(\mathbf{F}\_i^H \mathbf{F}\_{\rm BB}^{(n)} \mathbf{F}\_{\rm RF,i}\right) + \mathbf{F}\_{\rm RF,i}^H \mathbf{F}\_{\rm BB}^{(n)} \mathbf{F}\_{\rm BB}^{(n)H} \mathbf{F}\_{\rm RF,i} \tag{87}$$

$$s.t. \,\mathsf{C}\_1: \mathsf{F}\_{\mathsf{RF}} \in \mathbb{F}\_\prime \tag{88}$$

where *F H i* represents the *i*th row vector of the matrix *F*, and *F H RF*,*i* represents the *i*th row vector of the matrix *FRF*.

It should be noted that the third term *F H RF*,*i F* (*n*) *BB F* (*n*) *BB <sup>H</sup>FRF*,*<sup>i</sup>* in (86) is a convex function term, which needs further conversion. According to the first-order Taylor expansion, *F H RF*,*i F* (*n*) *BB F* (*n*) *BB <sup>H</sup>FRF*,*<sup>i</sup>* can be converted into:

$$\begin{aligned} \mathbf{F}\_{\mathrm{RF},i}^{H} \mathbf{F}\_{\mathrm{BB}}^{(n)} \mathbf{F}\_{\mathrm{BB}}^{(n)} \mathbf{F}\_{\mathrm{RF},i} &= \mathbf{F}\_{\mathrm{RF},i}^{(q)H} \mathbf{F}\_{\mathrm{BB}}^{(n)} \mathbf{F}\_{\mathrm{BB}}^{(n)} \mathbf{F}\_{\mathrm{RF},i}^{(q)} + 2\mathfrak{R} \left( \mathbf{F}\_{\mathrm{RF},i}^{(q)H} \mathbf{F}\_{\mathrm{BB}}^{(n)} \mathbf{F}\_{\mathrm{BB}}^{(n)} \mathbf{F}\_{\mathrm{BB}}^{(n)H} \left( \mathbf{F}\_{\mathrm{RF},i} - \mathbf{F}\_{\mathrm{RF},i}^{(q)} \right) \right) \\ &+ \left( \mathbf{F}\_{\mathrm{RF},i} - \mathbf{F}\_{\mathrm{RF},i}^{(q)} \right)^{H} \mathbf{F}\_{\mathrm{BB}}^{(n)} \mathbf{F}\_{\mathrm{BB}}^{(n)} \mathbf{F}\_{\mathrm{BB}}^{(n)} \left( \mathbf{F}\_{\mathrm{RF},i} - \mathbf{F}\_{\mathrm{RF},i}^{(q)} \right) \end{aligned} \tag{89}$$

where *F* (*q*) *RF*,*i* represents the estimated value of *FRF*,*<sup>i</sup>* at the *q*th iteration. According to the MM algorithm, the surrogate of *F H RF*,*i F* (*n*) *BB F* (*n*) *BB <sup>H</sup>FRF*,*<sup>i</sup>* can be expressed as follows:

$$\begin{split} \mathbf{F}\_{\textrm{RF},i}^{H} \mathbf{F}\_{\textrm{BB}}^{(n)} \mathbf{F}\_{\textrm{BB}}^{(n)} \mathbf{F}\_{\textrm{RF},i} \leq \mathbf{F}\_{\textrm{RF},i}^{(q)H} \mathbf{F}\_{\textrm{BB}}^{(n)} \mathbf{F}\_{\textrm{BB}}^{(n)} \mathbf{F}\_{\textrm{RF},i}^{(q)} + 2\mathfrak{R} \left( \mathbf{F}\_{\textrm{RF},i}^{(q)H} \mathbf{F}\_{\textrm{BB}}^{(n)} \mathbf{F}\_{\textrm{BB}}^{(n)} \mathbf{F}\_{\textrm{BB}}^{(n)H} \left( \mathbf{F}\_{\textrm{RF},i} - \mathbf{F}\_{\textrm{RF},i}^{(q)} \right) \right) \\ + \left( \mathbf{F}\_{\textrm{RF},i} - \mathbf{F}\_{\textrm{RF},i}^{(q)} \right)^{H} \mathbf{X}^{(n)} \left( \mathbf{F}\_{\textrm{RF},i} - \mathbf{F}\_{\textrm{RF},i}^{(q)} \right) \end{split} \tag{90}$$

where *X* (*n*) is a positive semidefinite matrix and satisfies the constraint *X* (*n*)*<sup>F</sup>* (*n*) *BB F* (*n*) *BB H*; here, we let *X* (*n*) <sup>=</sup> *<sup>λ</sup>*max *F* (*n*) *BB F* (*n*) *BB H <sup>I</sup>* and *<sup>λ</sup>*max *F* (*n*) *BB F* (*n*) *BB H* represents the maximum eigenvalue of the matrix *F* (*n*) *BB F* (*n*) *BB <sup>H</sup>*. In conclusion, (89) can be further expressed as:

$$\begin{split} \mathbf{F}\_{\mathrm{RF},i}^{H} \mathbf{F}\_{\mathrm{BB}}^{(n)} \mathbf{F}\_{\mathrm{BB}}^{(n)H} \mathbf{F}\_{\mathrm{RF},i} &\leq \lambda\_{\max} \Big( \mathbf{F}\_{\mathrm{BB}}^{(n)} \mathbf{F}\_{\mathrm{BB}}^{(n)H} \Big) \mathbf{F}\_{\mathrm{RF},i}^{H} \mathbf{F}\_{\mathrm{RF},i} + 2 \Re \Big( \mathbf{F}\_{\mathrm{RF},i}^{H} \Big( \mathbf{F}\_{\mathrm{BB}}^{(n)} \mathbf{F}\_{\mathrm{BB}}^{(n)H} - \lambda\_{\max} \Big( \mathbf{F}\_{\mathrm{BB}}^{(n)} \mathbf{F}\_{\mathrm{BB}}^{(n)H} \Big) \mathbf{I} \Big) \mathbf{F}\_{\mathrm{RF},i} \\ &\quad + \mathbf{F}\_{\mathrm{RF},i}^{(q)H} \Big( \lambda\_{\max} \Big( \mathbf{F}\_{\mathrm{BB}}^{(n)} \mathbf{F}\_{\mathrm{BB}}^{(n)H} \Big) \mathbf{I} - \mathbf{F}\_{\mathrm{BB}}^{(n)} \mathbf{F}\_{\mathrm{BB}}^{(n)H} \Big) \mathbf{F}\_{\mathrm{RF},i} \end{split} \tag{91}$$

According to (90), the surrogate function of the objective function in the problem *P*<sup>3</sup> can be expressed as:

$$\begin{aligned} P\_3: \min\_{\begin{subarray}{c} F\_{RF} \end{subarray}} & \sum\_{i=1}^{N} F\_i^H F\_i - 2\mathfrak{R} \Big( F\_i^H F\_{BB}^{(n)} F\_{RF,i} \Big) + F\_{RF,i}^H F\_{BB}^{(n)} F\_{BB}^{(n)H} F\_{RF,i} \Rightarrow \\ P\_4: \min\_{\begin{subarray}{c} F\_{RF} \end{subarray}} & \sum\_{i=1}^{N} F\_i^H F\_i - 2\mathfrak{R} \Big( F\_i^H F\_{BB}^{(n)} F\_{RF,i} \Big) + \lambda \Big( F\_{BB}^{(n)} F\_{BB}^{(n)H} \Big)\_{RF,i}^H \Big)\_{RF,i}^H \text{r}F\_i. \max \\ & \qquad + 2\mathfrak{R} (F\_{RF,i}^H F\_{BB}^{(n)} F\_{BB}^{(n)H}) - \\ & \lambda (F\_{BB}^{(n)} F\_{BB}^{(n)H})\_{max} (\lambda\_{RF,i}) (\binom{q}{RF,i}^H (\lambda \{F\_{BB}^{(n)} F\_{BB}^{(n)H}\}\_{BB}^{(n)} (\stackrel{n}{)} (\!)^{(n)})\_{\text{r}}^H))\_i \end{aligned} \tag{92}$$

It is worth noting that the first and third terms of the objective function in the problem *P*<sup>4</sup> are constant terms, and the last term is independent of the variable *F H RF*,*i* . After ignoring the above three items, the problem *P*<sup>4</sup> can be converted to the following projection problem:

$$P\_{\mathsf{S}} : \min\_{\mathsf{F}\_{RF}} \sum\_{i=1}^{N} \left\| \mathbf{F}\_{RF,i} - \mathbf{c}\_{i}^{(q)} \right\|\_{2'}^{2} \tag{93}$$

$$s.t. \,\mathsf{C}\_1: \mathsf{F}\_{\mathsf{RF}} \in \mathbb{F}\_\prime \tag{94}$$

where *c* (*q*) *<sup>i</sup>* = *F* (*n*) *BB F<sup>i</sup>* − *F* (*n*) *BB F* (*n*) *BB <sup>H</sup>* <sup>−</sup> *<sup>λ</sup>*max *F* (*n*) *BB F* (*n*) *BB H I F* (*q*)*H RF*,*i* .

Therefore, the following closed form solution can be obtained for the problem *P*5, which is reflected in line three of the **inner algorithm** in Algorithm 3:

$$F\_{RF,i} = e^{\text{jarg}(c\_i^{(q)})} \, \_\prime \forall i \,, \tag{95}$$

$$F\_{RF} = e^{-j \text{arg}(\mathbb{C}^{(q)T})},\tag{96}$$

where *C* (*q*) = *F* (*n*) *BB F <sup>H</sup>* <sup>−</sup> *F* (*n*) *BB F* (*n*) *BB <sup>H</sup>* <sup>−</sup> *<sup>λ</sup>*max *F* (*n*) *BB F* (*n*) *BB H I F* (*q*)*H RF* , which is reflected in line two of the **inner algorithm** in Algorithm 3.

#### *7.4. Solution of The Digital Beamforming Matrix*

Based on the analog beamforming matrix *FRF* obtained in the previous section, the solution problem of the digital beamforming matrix *FBB* can be modeled as follows:

$$P\_6: \min\_{F\_{BB}} \left\| F - F\_{RF}^{(n)} F\_{BB} \right\|\_{F'}^2 \tag{97}$$

$$\text{s.t.} \left\| \mathbf{F}\_{\text{RF}}^{(n)} \mathbf{F}\_{\text{BB}} \right\|\_{F}^{2} = \text{P}\_{\text{T}} \tag{98}$$

where *F* (*n*) *RF* represents the estimated value of the analog beamforming matrix *FRF* at the *n*th iteration. Due to the quadratic form and the convex equality constraint in the problem *P*6, the problem *P*<sup>6</sup> is a nonconvex QCQP form problem.

One of the ways to solve the problem *P*<sup>6</sup> is to relax the equality constraint into the inequality constraint, and then convert the problem *P*<sup>6</sup> into a convex minimization problem, which can be solved with the CVX toolbox, but the complexity of this method is high. To this end, based on the fact that the hybrid beamforming matrix *F* satisfies the power constraint, i.e., k*F*k 2 *<sup>F</sup>* = *PT*, the following closed form solution can be obtained for the problem *P*6, which is reflected in line two of the **main algorithm** in Algorithm 3:

$$\mathbf{F}\_{BB} = \left(\mathbf{F}\_{RF}^{(n)H} \mathbf{F}\_{RF}^{(n)}\right)^{-1} \mathbf{F}\_{RF}^{(n)H} \mathbf{F}\_{\prime} \tag{99}$$

In conclusion, the MM-AltOpt algorithm for solving *FBB* and *FRF* can be described as the main algorithm and the inner algorithm, as follows:

**Algorithm 3:** Design algorithm of the digital beamforming matrix and the analog beamforming matrix.

#### **Main algorithm: MM-AltOpt algorithm.**

**Input**: Hybrid beamforming matrix *F*, initial: *F* (*n*=0) *RF* ∈ F, iteration index *n* = 0, threshold *ε*<sup>3</sup> = 10−<sup>3</sup> , the solution of the objective function of the problem *P*<sup>1</sup> in the *n*th iteration is *δ* (*n*) . 1. while *δ* (*n*) <sup>−</sup> *<sup>δ</sup>* (*n*−1) <sup>≥</sup> *<sup>ε</sup>*<sup>3</sup>


5. end

**Output**: *FRF*,*FBB*, normalize *FBB* = √ *PT* k*FRF <sup>F</sup>BB*k*<sup>F</sup> FBB*.

**Inner algorithm: Algorithm for solving the analog beamforming matrix.**

**Input**: Hybrid beamforming matrix *F*, *F* (*n*) *BB* , *F* (*q*=0) *RF* <sup>∈</sup> <sup>F</sup>, iteration index *<sup>q</sup>* <sup>=</sup> 0, threshold *<sup>ε</sup>*<sup>4</sup> <sup>=</sup> <sup>10</sup>−<sup>3</sup> , the solution of the objective function of the problem *P*<sup>5</sup> in the *q*th iteration is *δ* (*q*) 1 .

```
1. while
            
            
            
             δ
              (n) − δ
                       (n−1)
                              
                              
                               ≥ ε4
```

$$\text{2. Calculate } \overset{\cdot}{\mathbb{C}}^{(q)} = \overset{\cdot}{\mathbf{F}\_{BB}^{(n)}} \overset{\cdot}{\mathbf{F}^{H}} - \left(\overset{\cdot}{\mathbf{F}\_{BB}^{(n)}} \overset{\cdot}{\mathbf{F}\_{BB}^{(n)}} - \lambda\_{\text{max}} \left(\mathbf{F}\_{BB}^{(n)} \overset{\cdot}{\mathbf{F}\_{BB}^{(n)}} \overset{\cdot}{\mathbf{F}\_{BB}^{(n)}}\right) \mathbf{I}\right) \overset{(q)H}{\underset{RF}{\mathbf{F}}^{(q)}} \overset{\cdot}{\mathbf{F}}$$

.

3. Calculate *FRF* = *e* −*j*arg(*C* (*q*)*T* ) 4. Set *q* = *q* + 1.

5. end

**Output**: *F* (*n*) *RF* .

#### **8. Results and Discussion**

In this section, we evaluate the performance of the proposed joint user scheduling and hybrid beamforming design algorithm by numerical simulations. In the numerical simulations, we set the number of multicast groups to *L* = 7, which cover 150 active users, and set the SINR constraint threshold of each multicast group to SINR<sup>0</sup> = 1. To facilitate analysis, we assume that the CSI errors of different multicast groups are the same, which

are expressed as *σ* 2 *f rsd d*,*k*,*l* = *σ* 2 *f rsd d* and *σ* 2 *φτk*,*<sup>l</sup>* = *σ* 2 *φτ* . The value of *P*<sup>0</sup> can be calculated by [31]. In addition, the system parameters used in the numerical simulations are shown in Table 1.

**Table 1.** Simulation parameters.


Figure 7 shows the convergence trajectory of the EE of the massive MIMO LEO sat-

Figure 7 shows the convergence trajectory of the EE of the massive MIMO LEO satellite multigroup multicast communication system, versus the number of iterations for different CSI errors, different numbers of preselected users and different scheduling algorithms. In this simulation, two groups of channel errors are set according to Refs. [23,29], i.e., *σ* 2 *f rsd d* = 25, *σ* 2 *φτ* = 10 and *σ* 2 *f rsd d* = 20, *σ* 2 *φτ* = 5. In addition, we set *U<sup>s</sup>* = 2 and *P<sup>T</sup>* = 50 W. Meanwhile, we set two different numbers of the preselected users, i.e., *Ul*,*p*/*U<sup>s</sup>* = 2, *Ul*,*<sup>p</sup>* = 4 and *Ul*,*p*/*U<sup>s</sup>* = 3, *Ul*,*<sup>p</sup>* = 6. It can be seen that the proposed robust algorithm has higher performance gain than the traditional nonrobust algorithm, which shows the effectiveness of the robust algorithm. Meanwhile, it can be seen that when *σ* 2 *f rsd d* = 25, *σ* 2 *φτ* = 10 and *σ* 2 *f rsd d* = 20, and *σ* 2 *φτ* = 5, the EE performance gain of the proposed robust algorithm is improved by 9.8% and 6.7%, respectively, compared with the traditional nonrobust algorithm. The system EE of the proposed joint user scheduling and hybrid beamforming design algorithm is higher than that of the decoupling design algorithm, and the more preselected users, the higher performance improvement. ellite multigroup multicast communication system, versus the number of iterations for different CSI errors, different numbers of preselected users and different scheduling algorithms. In this simulation, two groups of channel errors are set according to Refs. [23,29], i.e., 2 *rsd* 25 *d f* = , 2 10 = and 2 *rsd* 20 *d f* = , 2 5 = . In addition, we set 2 *U<sup>s</sup>* = and *PT* = 50W . Meanwhile, we set two different numbers of the preselected users, i.e., , , / 2, 4 *U U U l p s l p* = = and , , / 3, 6 *U U U l p s l p* = = . It can be seen that the proposed robust algorithm has higher performance gain than the traditional nonrobust algorithm, which shows the effectiveness of the robust algorithm. Meanwhile, it can be seen that when 2 *rsd* 25 *d f* = , 2 10 = and 2 *rsd* 20 *d f* = , and 2 5 = , the EE performance gain of the proposed robust algorithm is improved by 9.8% and 6.7%, respectively, compared with the traditional nonrobust algorithm. The system EE of the proposed joint user scheduling and hybrid beamforming design algorithm is higher than that of the decoupling design algorithm, and the more preselected users, the higher performance improvement.

**Figure 7.** Convergence trajectory of system EE relative to different CSI errors, different number of preselected users and different scheduling algorithms. preselected users and different scheduling algorithms.

Figure 8 compares the EE performance of the proposed algorithm and the traditional algorithm under different system parameters, versus different transmission power thresholds *PT* . It can be seen that with the increase in transmission power, the EE performance shows a trend of first rising and then falling. The reason is that the growth rate of the system rate is lower than that of the power consumption. Meanwhile, we can see that the 28.41% and 45.19% higher than that of the traditional decoupling design algorithm. **Figure 7.** Convergence trajectory of system EE relative to different CSI errors, different number of Figure 8 compares the EE performance of the proposed algorithm and the traditional algorithm under different system parameters, versus different transmission power thresholds *PT*. It can be seen that with the increase in transmission power, the EE performance shows a trend of first rising and then falling. The reason is that the growth rate of the system rate is lower than that of the power consumption. Meanwhile, we can see that the

EE of the joint user scheduling and hybrid beamforming design algorithm is higher than that of the decoupling design algorithm versus different transmission power thresholds

> and *PT* =15W

, the EE performance gain of the proposed joint design algorithm is

, when

, , / 2, 4 *U U U l p s l p* = =

2 10 

=

2 *rsd* 25 *d f* = ,

*PT*

and

, , / 3, 6 *U U U l p s l p* = =

EE of the joint user scheduling and hybrid beamforming design algorithm is higher than that of the decoupling design algorithm versus different transmission power thresholds *PT*. Under the conditions of *σ* 2 *f rsd d* = 25, *σ* 2 *φτ* = 10 and *P<sup>T</sup>* = 15 W, when *Ul*,*p*/*U<sup>s</sup>* = 2, *Ul*,*<sup>p</sup>* = 4 and *Ul*,*p*/*U<sup>s</sup>* = 3, *Ul*,*<sup>p</sup>* = 6, the EE performance gain of the proposed joint design algorithm is 28.41% and 45.19% higher than that of the traditional decoupling design algorithm. *Sensors* **2022**, *22*, x FOR PEER REVIEW 27 of 32

**Figure 8.** Comparison of system EE of different algorithms with different transmission power thresholds *PT* . **Figure 8.** Comparison of system EE of different algorithms with different transmission power thresholds *PT*.

Figure 9 indicates the change trend of the SE of the massive MIMO LEO satellite multigroup multicast communication system versus different transmission power thresholds *PT* . It can be seen that the SE increases with the increase in the transmission power. Meanwhile, we can see that the SE of the joint user scheduling and hybrid beamforming design algorithm is higher than that of the decoupling design algorithm. In addition, the more preselected users, the higher the system SE. Compared with the traditional algorithm, the proposed robust joint design algorithm can obtain higher system SE at the same transmission power, and thus can improve the system EE. Under the conditions of 2 *rsd* 25 *d f* = , 2 10 = and *PT* = 30W , when , , / 2, 4 *U U U l p s l p* = = and , , / 3, 6 *U U U l p s l p* = = , with the improvement of system SE performance, the EE performance gain of the proposed joint design algorithm is 26.16% and 37.85% higher than that of the traditional decoupling design al-Figure 9 indicates the change trend of the SE of the massive MIMO LEO satellite multigroup multicast communication system versus different transmission power thresholds *PT*. It can be seen that the SE increases with the increase in the transmission power. Meanwhile, we can see that the SE of the joint user scheduling and hybrid beamforming design algorithm is higher than that of the decoupling design algorithm. In addition, the more preselected users, the higher the system SE. Compared with the traditional algorithm, the proposed robust joint design algorithm can obtain higher system SE at the same transmission power, and thus can improve the system EE. Under the conditions of *σ* 2 *f rsd d* = 25, *σ* 2 *φτ* = 10 and *P<sup>T</sup>* = 30 W, when *Ul*,*p*/*U<sup>s</sup>* = 2, *Ul*,*<sup>p</sup>* = 4 and *Ul*,*p*/*U<sup>s</sup>* = 3, *Ul*,*<sup>p</sup>* = 6, with the improvement of system SE performance, the EE performance gain of the proposed joint design algorithm is 26.16% and 37.85% higher than that of the traditional decoupling design algorithm.

gorithm. Figure 10 shows the SE comparison of different multicast groups. It can be seen that the SE of each multicast group of the proposed robust algorithm is higher than that of the nonrobust algorithm. Meanwhile, with the increase in the number of preselected users, the diversity of users increases, and the performance of the proposed joint user scheduling and hybrid beamforming design algorithm also improves. This is because with the increase in the number of preselected users, the range of users that can be scheduled and selected increases. By scheduling different users in each multicast group, the SE can be further improved. Meanwhile, with the improvement of the system SE, the system EE performance gain also increases, which verifies the effectiveness of the joint user scheduling and hybrid beamforming design algorithm.

**Figure 9.** Change trajectory of system SE versus different transmission power thresholds

*PT* . thresholds

*PT*

2 10 

gorithm.

= *PT* .

and *PT* = 30W , when

**Figure 9.** Change trajectory of system SE versus different transmission power thresholds *PT* **Figure 9.** Change trajectory of system SE versus different transmission power thresholds *P* . *T*. ing and hybrid beamforming design algorithm.

**Figure 8.** Comparison of system EE of different algorithms with different transmission power

Figure 9 indicates the change trend of the SE of the massive MIMO LEO satellite multigroup multicast communication system versus different transmission power thresholds

. It can be seen that the SE increases with the increase in the transmission power. Meanwhile, we can see that the SE of the joint user scheduling and hybrid beamforming design algorithm is higher than that of the decoupling design algorithm. In addition, the more preselected users, the higher the system SE. Compared with the traditional algorithm, the proposed robust joint design algorithm can obtain higher system SE at the same transmis-

ment of system SE performance, the EE performance gain of the proposed joint design algorithm is 26.16% and 37.85% higher than that of the traditional decoupling design al-

and

, , / 3, 6 *U U U l p s l p* = =

2 *rsd* 25 *d f* = ,

, with the improve-

sion power, and thus can improve the system EE. Under the conditions of

, , / 2, 4 *U U U l p s l p* = =

Figure 11 shows the change trajectory of the system EE and SE versus the different *Us* . It can be seen that with the increase in *Us* , the system EE and SE show a downward trend. This is because the communication rate of each multicast group is constrained by the user with the worst SINR in the multicast group. With the increase in *Us* , if the users' channel vectors in the multicast group remain collinear, the EE and SE would remain unchanged. However, according to the rules of the user preselection and scheduling, with the increase in *Us* , the collinearity among users in the multicast group would decrease, which can result in the increase in interference and the decrease in the worst SINR in each multicast group. In other words, with the decrease in *Us* , the users**'** SINR will be improved. Therefore, with the improvement of SINR, the system EE performance gain also increases, as shown in Figure 11a. Under the condition of / 2, 4 *U U U l p s l p* = = , when 2 *U<sup>s</sup>* = Figure 11 shows the change trajectory of the system EE and SE versus the different *U<sup>s</sup>* . It can be seen that with the increase in *U<sup>s</sup>* , the system EE and SE show a downward trend. This is because the communication rate of each multicast group is constrained by the user with the worst SINR in the multicast group. With the increase in *U<sup>s</sup>* , if the users' channel vectors in the multicast group remain collinear, the EE and SE would remain unchanged. However, according to the rules of the user preselection and scheduling, with the increase in *U<sup>s</sup>* , the collinearity among users in the multicast group would decrease, which can result in the increase in interference and the decrease in the worst SINR in each multicast group. In other words, with the decrease in *U<sup>s</sup>* , the users' SINR will be improved. Therefore, with the improvement of SINR, the system EE performance gain also increases, as shown in Figure 11a. Under the condition of *Ul*,*p*/*U<sup>s</sup>* = 2, *Ul*,*<sup>p</sup>* = 4, when *U<sup>s</sup>* = 2, the EE performance gain is 21.05% higher than when *U<sup>s</sup>* = 4.

, ,

4 *U<sup>s</sup>* = .

, the EE performance gain is 21.05% higher than when

*Sensors* **2022**, *22*, x FOR PEER REVIEW 29 of 32

Figure 12 shows the performance comparison of different algorithms for solving *FBB* and *FRF* . In this simulation, we set three comparison algorithms, i.e., the optimal design algorithm, the alternating minimization algorithm based on the phase extraction (PE-Altmin) algorithm and the orthogonal matching pursuit (OMP) algorithm. The optimal design algorithm refers to the numerical simulation result of the hybrid beamforming matrix *F* . It can be seen that the system performance of the MM-AltOpt algorithm is slightly lower than that of the optimal design algorithm. Meanwhile, in Figure 12a, when *PT* =10W , we can see that the system EE performance gain of the proposed MM-AltOpt algorithm is improved by about 2% and 5%, respectively, compared with the PE-Altmin algorithm and the OMP algorithm. In addition, from the perspective of algorithm complexity, the complexities of the MM-AltOpt algorithm, PE-Altmin algorithm and OMP algorithm are ( ( )) 3 2 <sup>2</sup> *O I N I NL NL MM* + + *Inner* , ( ( )) 3 3 *O I N L NL PE* + + and ( ( )) 4 2 2 2 3 <sup>2</sup> *O I L N L N L L OMP* + + + , respectively, where *MM I* , *Inner I* , *PE I* and *OMP I* are the number of iterations of the corresponding algorithm. In conclusion, we can see that the complexity of the proposed MM-AltOpt algorithm is close to that of the other two algorithms, however, the system EE performance is higher, which can verify the effectiveness of the proposed algorithm. Figure 12 shows the performance comparison of different algorithms for solving *FBB* and *FRF*. In this simulation, we set three comparison algorithms, i.e., the optimal design algorithm, the alternating minimization algorithm based on the phase extraction (PE-Altmin) algorithm and the orthogonal matching pursuit (OMP) algorithm. The optimal design algorithm refers to the numerical simulation result of the hybrid beamforming matrix *F*. It can be seen that the system performance of the MM-AltOpt algorithm is slightly lower than that of the optimal design algorithm. Meanwhile, in Figure 12a, when *P<sup>T</sup>* = 10 W, we can see that the system EE performance gain of the proposed MM-AltOpt algorithm is improved by about 2% and 5%, respectively, compared with the PE-Altmin algorithm and the OMP algorithm. In addition, from the perspective of algorithm complexity, the complexities of the MM-AltOpt algorithm, PE-Altmin algorithm and OMP algorithm are *O IMM N*<sup>3</sup> + *IInner*2*NL*<sup>2</sup> + *NL*, *O IPE N*<sup>3</sup> + *L* <sup>3</sup> + *NL* and *O IOMP L* <sup>4</sup>*N* + *L* <sup>2</sup> + *N*2*L* <sup>2</sup> + 2*L* 3 , respectively, where *IMM*, *IInner*, *IPE* and *IOMP* are the number of iterations of the corresponding algorithm. In conclusion, we can see that the complexity of the proposed MM-AltOpt algorithm is close to that of the other two algorithms, however, the system EE performance is higher, which can verify the effectiveness of the proposed algorithm. Figure 12 shows the performance comparison of different algorithms for solving *FBB* and *FRF* . In this simulation, we set three comparison algorithms, i.e., the optimal design algorithm, the alternating minimization algorithm based on the phase extraction (PE-Altmin) algorithm and the orthogonal matching pursuit (OMP) algorithm. The optimal design algorithm refers to the numerical simulation result of the hybrid beamforming matrix *F* . It can be seen that the system performance of the MM-AltOpt algorithm is slightly lower than that of the optimal design algorithm. Meanwhile, in Figure 12a, when *PT* =10W , we can see that the system EE performance gain of the proposed MM-AltOpt algorithm is improved by about 2% and 5%, respectively, compared with the PE-Altmin algorithm and the OMP algorithm. In addition, from the perspective of algorithm complexity, the complexities of the MM-AltOpt algorithm, PE-Altmin algorithm and OMP algorithm are ( ( )) 3 2 <sup>2</sup> *O I N I NL NL MM* + + *Inner* , ( ( )) 3 3 *O I N L NL PE* + + and ( ( )) 4 2 2 2 3 <sup>2</sup> *O I L N L N L L OMP* + + + , respectively, where *MM I* , *Inner I* , *PE I* and *OMP I* are the number of iterations of the corresponding algorithm. In conclusion, we can see that the complexity of the proposed MM-AltOpt algorithm is close to that of the other two algorithms, however, the system EE performance is higher, which can verify the effectiveness of the proposed algorithm.

**Figure 12.** Comparison of system performance of different algorithms for solving *FBB* and *FRF*. (**a**) Comparison of system EE of different algorithms for solving *FBB* and *FRF*; (**b**) comparison of system SE of different algorithms for solving *FBB* and *FRF*.

#### **9. Conclusions**

In this paper, we investigated the robust joint user scheduling and hybrid beamforming design scheme for the massive MIMO LEO satellite multigroup multicast communication system. The scheme design considered the limited transmission power of the LEO satellite and the requirement of QoS and analyzed the influence of residual Doppler shift and phase disturbance on CSI errors. On this basis, taking the system EE as the optimization objective, we focused on the robust joint user scheduling and hybrid beamforming design. To reduce the complexity of the algorithm, we proposed the user preselection step, which can significantly reduce the system complexity while ensuring the system performance. For the nonconvex problem of the objective function, we adopted the CCCP framework after transforming the optimization problem into the DC programming problem. For the rank-one constraint, we proposed the penalty iterative algorithm. Finally, to obtain the digital and analog beamforming matrices, we adopted the MM-AltOpt algorithm.

Numerical results indicated that the proposed algorithm can effectively improve the system EE. The EE performance gain of the proposed robust algorithm was improved by nearly 10% compared with the traditional nonrobust algorithm. Meanwhile, the EE performance gain of the proposed joint user scheduling and hybrid beamforming design algorithm was improved by nearly 40% compared with the traditional decoupling design algorithm. In conclusion, the robust joint user scheduling and hybrid beamforming design algorithm proposed in this paper can significantly improve the system EE performance.

**Author Contributions:** Conceptualization, Y.L. and C.L.; methodology, Y.L. and J.L.; software, Y.L. and L.F.; validation, Y.L., C.L. and J.L.; formal analysis, Y.L. and C.L.; investigation, Y.L.; writing original draft preparation, Y.L. and L.F.; writing—review and editing, C.L. and J.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by National Natural Science Foundation of China, grant number 62001516.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available in the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**



#### **References**


**Mayada Osama 1,2,\*, Salwa El Ramly <sup>2</sup> and Bassant Abdelhamid <sup>2</sup>**


**Abstract:** The dense deployment of small cells (SCs) in the 5G heterogeneous networks (HetNets) fulfills the demand for vast connectivity and larger data rates. Unfortunately, the power efficiency (PE) of the network is reduced because of the elevated power consumption of the densely deployed SCs and the interference that arise between them. An approach to ameliorate the PE is proposed by switching off the redundant SCs using machine learning (ML) techniques while sustaining the quality of service (QoS) for each user. In this paper, a linearly increasing inertia weight–binary particle swarm optimization (IW-BPSO) algorithm for SC on/off switching is proposed to minimize the power consumption of the network. Moreover, a soft frequency reuse (SFR) algorithm is proposed using classification trees (CTs) to alleviate the interference and elevate the system throughput. The results show that the proposed algorithms outperform the other conventional algorithms, as they reduce the power consumption of the network and the interference among the SCs, ameliorating the total throughput and the PE of the system.

**Keywords:** 5G HetNets; BPSO; classification trees (CTs); soft frequency reuse (SFR); small cells (SCs)

## **1. Introduction**

Recently, the exponential growth of the numerous wireless devices and the datahungry applications have earned huge significance. This required imperious expansion of the 5G network to support the forthcoming 5G use cases, such as video live-streaming, conferencing, online gaming, etc. [1,2]. Moreover, the 5G cellular network is planned to elevate the capacity 1000 times and the spectrum efficiency by 5–15 times with respect to 4G [3,4]. This can be achieved by utilizing heterogeneous networks (HetNets), which can enhance the system data rates and the quality of service (QoS) of the users as the small cells (SCs) are deployed within the macro cells (MCs) coverage area. Furthermore, SCs offer the benefit of providing service to previously uncovered regions and in the network regions demanding larger capacity [2,3,5–7]. Figure 1 shows a general representation of the HetNet scenario with the MCs underlaid by the densely deployed SCs.

It is foreseen that the massive growth in the SC deployment will be continued in the coming years [8,9], leading to various challenges such as the interference among the SCs [10], their elevated power consumption [11,12], and the elevated operating expenses [9]. Thus, it is crucial to face these challenges to ameliorate the performance of 5G HetNets.

**Citation:** Osama, M.; El Ramly, S.; Abdelhamid, B. Binary PSO with Classification Trees Algorithm for Enhancing Power Efficiency in 5G Networks. *Sensors* **2022**, *22*, 8570. https://doi.org/10.3390/s22218570

Academic Editor: Josip Lorincz

Received: 5 October 2022 Accepted: 4 November 2022 Published: 7 November 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

**Figure 1.** General representation of HetNet scenario with densely deployed small cells. **Figure 1.** General representation of HetNet scenario with densely deployed small cells.

It is foreseen that the massive growth in the SC deployment will be continued in the coming years [8, 9], leading to various challenges such as the interference among the SCs [10], their elevated power consumption [11, 12], and the elevated operating expenses [9]. Thus, it is crucial to face these challenges to ameliorate the performance of 5G HetNets. The main objective of our paper is to propose a new approach for the irregular nature of the 5G HetNets that merge the usage of both binary particle swarm optimization (BPSO) algorithm with linear increasing inertia weight (IW) and soft frequency reuse (SFR) to maximize the power efficiency (PE) of the SCs and minimize the number of active SCs while guaranteeing the QoS for the UEs. In SFR, every SC is split into center and edge regions where one of the unutilized sub-bands by the edge regions of the adjoining SCs is allocated to the edge region of the SC, while the remaining sub-bands are utilized in the center region of the SC with reduced transmission power to diminish the interference to the adjoining SCs. Unlike the prior works [9, 13, 14] that will be discussed later in Section 2, the proposed algorithm utilizes the linear increasing IW-BPSO algorithm for selecting the SCs to be switched on/off, and we propose the SFR utilizes classification trees (CTs) to The main objective of our paper is to propose a new approach for the irregular nature of the 5G HetNets that merge the usage of both binary particle swarm optimization (BPSO) algorithm with linear increasing inertia weight (IW) and soft frequency reuse (SFR) to maximize the power efficiency (PE) of the SCs and minimize the number of active SCs while guaranteeing the QoS for the UEs. In SFR, every SC is split into center and edge regions where one of the unutilized sub-bands by the edge regions of the adjoining SCs is allocated to the edge region of the SC, while the remaining sub-bands are utilized in the center region of the SC with reduced transmission power to diminish the interference to the adjoining SCs. Unlike the prior works [9,13,14] that will be discussed later in Section 2, the proposed algorithm utilizes the linear increasing IW-BPSO algorithm for selecting the SCs to be switched on/off, and we propose the SFR utilizes classification trees (CTs) to enhance the network PE, taking into consideration the irregular nature of the 5G HetNets. The on/off switching of SCs using BPSO algorithm is carried out first; then the SFR is applied for sub-band allocation to minimize the number of operations required for allocating the sub-bands to the SCs.

enhance the network PE, taking into consideration the irregular nature of the 5G HetNets. The main contributions of this paper can be summarized as follow:


vergence of the BPSO algorithm. 2. Propose a novel frequency allocation algorithm for SFR based on the CTs as it is simple and accurate machine learning (ML) technique to mitigate the interference among Results demonstrate that the proposed algorithms have superior performance over the conventional algorithms (always on, random 10%, and BPSO only), as it has higher total system throughput and PE, and lower system power consumption and outage probability.

the irregularly shaped SCs. Results demonstrate that the proposed algorithms have superior performance over The remainder of the paper is organized as follows: the literature review is presented in Section 2. Section 3 demonstrates the system model, and Section 4 explains the proposed

the conventional algorithms (always on, random 10%, and BPSO only), as it has higher

algorithms. Then, Section 5 shows the simulation results. Eventually, Section 6 concludes the paper.

#### **2. Literature Review**

Recently, immense research has been carried out to ameliorate the PE of the mod-ern communication networks, such as satellite and terrestrial networks [15,16], massive MIMO systems [17] and SCs networks [18–23].

SC on/off switching is an auspicious approach to minimize the power consumption and enhance the PE of the system [18–23]. The authors in [18,19] studied the elevated the energy consumption of WLAN. They proposed solving the problem by on/off switching and power adjustment of the access stations. The authors in [20] proposed a load-aware strategy, where SCs in HetNets are switched to sleep mode according to their load level. In [21], every SC independently switches off upon the decrease in the number of user equipments (UEs) and activates using one of three approaches. The first approach is the sleeping SC keeps sensing the interference plus noise levels and switches on when a new UE is sensed in its coverage area. The second approach is that the MC sends a wake-up request to all SCs upon increasing the number of UEs associated with the MC, then switches off the SCs with no UEs later. In the third approach, a time advance indicator is sent to the MC by the UEs and the SCs and is utilized by the MC to determine the nearest SC to the UE to switch on. In [22], pre-sleeping SCs at the same zone create a sleeping cluster. Then, the SCs in the sleeping cluster are randomly selected to be switched off leaving only one active SC to guarantee the coverage. On the other hand, the authors in [23] proposed switching off the SCs and handing over the UEs to the MC. However, the unplanned SC off switching may increase the unnecessary handovers and underutilize the SCs. Thus, novel SC on/off switching techniques are required to enhance the performance of the 5G HetNets and to reduce the power consumption of the system.

PSO is a prevalent meta-heuristic algorithm used in solving optimization problems [9]; thus, the authors in [9,24,25] utilized the PSO algorithm for SC switching to enhance the performance of the system. The authors in [9] proposed an efficient cell modeling (ECM) algorithm to set up the connection initially between the UEs and the SCs by selecting the strongest received signals. Then, the BPSO algorithm is utilized to turn off the excessive SCs. The authors in [24] proposed first utilization of BPSO algorithm to choose the MCs' optimum locations not only to achieve minimum overlap but also to guarantee a reasonable coverage for the UEs. Then, a multi-stage PSO (MS-PSO) algorithm consisting of two interactive loops are utilized. The outer loop is utilized to switch the SCs (on or off), while the inner loop is utilized to optimize the active power of the SC and to elevate the power of the SC if the data rate rises. On the other hand, a combined optimal frequency and power allocation (COFPA) scheme is proposed in [25]. First, using the BPSO algorithm, the MCs are switched on and off until the lead interferer is abolished with minimum cost function and with a reasonable coverage to the UEs. Then, the MS-PSO algorithm is utilized to control the SC switching to mitigate the interference and minimize the power consumption. However, the convergence of the PSO can be enhanced by adjusting its IW, leading to improve the system performance.

IW has a significant role in the process of offering a trade-off among diversification and intensification skills of PSO algorithm. Reducing the IW facilitates exploring the search space (global search), although raising the IW aids exploiting the search space (local search) [26] to find the solution (particle). Numerous approaches are presented to adjust the IW such as the constant IW [9,27] and the random IW [28]. However, the constant IW approach can fail to balance exploration and exploitation because of the lack of adjustment of IW [29,30]. On the other hand, the authors in [31] propose a linearly decreasing IW technique where the IW is initialized at a larger value; then it is linearly reduced to a smaller value. However, in this technique, the tendency of the particles to local search is constantly increasing. The authors in [32,33] demonstrated that the increase in IW surpasses the decrease in IW for PSO on all their tested benchmarks. Thus, the linearly increasing IW is adopted in our paper.

Various interference mitigation techniques for the modern networks are presented in the literature such as advanced multiple access techniques [34,35] and frequency reuse (FR) [13,36–39]. FR is an auspicious approach aiming to mitigate the interference in the modern HetNets [37]. SFR is presented to minimize the interference in HetNets [13,36,38,39]. The authors in [36] present SFR in HetNets, where the MC organizes the assignment of the resource plans to the SCs, and the SCs choose the resource plan. However, this algorithm cannot be implemented in the absence of the MC. A multi-level SFR (MSFR) for HetNets is demonstrated in [38], where every cell is split into three zones (central, intermediate, and edge), utilizing various frequency segments and transmission power levels. The authors in [13] proposed a novel SFR algorithm to mitigate the interference by splitting the SC to the center and edge regions. Moreover, the on/off switching of the SCs depends on their interference contribution rate (ICR) values. The authors in [39] presented an MSFR scheme where every MC is split into various circular areas, and different spectrum and power are allocated to every area. Unfortunately, in this scheme, SFR is not applied to the SCs but only the MCs. New SFR approaches are essential to alleviate the interference and improve the network throughput in 5G HetNets.

Moreover, the performance of 5G HetNets can be enhanced using ML techniques [40]. Decision trees (DTs) is one of the propitious supervised learning (SL) techniques, where every training example must be fed with its label to train a learning model, then utilize this model to predict the output for any new data. The authors in [41] proposed using supervised ML to enhance both the classification of services and the distribution of network resources in 5G networks. Moreover, the results reveal that DTs and random forests are the best approaches. DTs are split to classification trees (CTs) and regression trees [40]. The authors in [42] studied several anomaly detection techniques in 5G traffic. The performance of these techniques is analyzed based on multiple factors such as the probability of identifying anomalies and the probability of detecting a false positive. The results demonstrated that the CTs technique outperforms the other techniques. The problem of monitoring and predicting the quality of experience of cellular networks is studied in [43]. The authors compared various SL techniques and trained them utilizing training data based on traffic measurements of the UEs from a field trial study. The CTs are the ultimately chosen model because of their superb prediction accuracy and their prediction speed. In our paper, CT is utilized for the first time to allocate the frequencies in the SFR 5G network.

Finally, to match real-life scenarios, the irregular nature of the 5G HetNets was modeled utilizing Voronoi cells [44–47]. Because of the immense deployment of SCs, Voronoi cells are considered more practical than traditional hexagonal grids [48].

#### **3. System Model**

Consider 5G HetNet with densely deployed Voronoi SCs, where the SCs are deployed within the MCs coverage area. In this scenario, the MCs and the SCs utilize different frequency bands, mitigating the cross-tier interference between them. The MCs stay active to maintain the coverage of the network when the SCs are turned off. On the other hand, the SCs are either active (on) or asleep (off). When the SC is turned off, regular discovery signals are sent by the SC to be detected by any potential user. Each UE reports its channel state information and its reference signal received power (RSRP) to its SC. The SC sub-band allocation and the SC switching is organized by a main controller, or the MC if the main controller is absent, to collect the data from the SCs, allocate the sub-bands to the SCs, and determine the on/off switching decisions of the SCs. Furthermore, it is assumed that all UEs in every SC are located inside the coverage area of the SC. In case of the existence of any coverage holes in the SCs, relay nodes can be utilized to cover these holes [49]. However, this is not considered in our paper. Nomenclature lists the described symbols utilized in this paper.

An SC on/off indicator *γ* is defined, where *γ<sup>m</sup>* = 1 when the SC *m* is active; otherwise, *γ<sup>m</sup>* = 0. The number of SCs is denoted by "*M*" and the number of UEs in SC *m* is denoted by "*Um*". The UE association indicator *ϕm*,*<sup>k</sup>* = 1 if UE *k* is associated with SC *m*; otherwise, *ϕm*,*<sup>k</sup>* = 0. The signal to interference and noise ratio (SINR) of UE *k* in SC *m* can be calculated as [14]:

$$SINR\_{k,m} = \frac{\gamma\_m P\_m G\_{m,k}}{\sum\_{n \neq m,\ n \in M} \gamma\_n P\_n G\_{n,k} + N} \tag{1}$$

where *P<sup>m</sup>* and *P<sup>n</sup>* are the transmission powers of the serving SC *m* and the interfering SC *n*, respectively. The channel gain between UE *k* and serving SC *m* is *Gm*,*<sup>k</sup>* = *dm*,*<sup>k</sup>* −*α* , where *dm*,*<sup>k</sup>* is the distance between SC *m* and UE *k* and α is the path loss exponent [50]; the channel gain between UE *k* and interfering SC *n* is *Gn*,*<sup>k</sup>* and *N* is the noise power. The data rate *Rk*,*<sup>m</sup>* of UE *k* in SC *m* is also calculated by Shannon's formula as [13]:

$$R\_{k,m} = B\_k \, \log\_2(1 + SINR\_{k,m}) \tag{2}$$

where *B<sup>k</sup>* = *BRBLk*,*<sup>m</sup>* is the bandwidth allocated to UE *k*, while the resource block (RB) bandwidth is *BRB* and *Lk*,*<sup>m</sup>* is the number of requisite RBs for UE *k* in SC *m* to achieve the minimum data rate [14].

While the total throughput of the system is given by [9]:

$$\mathcal{C}\_{\text{sys}} = \sum\_{m=1}^{M} \gamma\_m \sum\_{k=1}^{U\_m} \mathcal{R}\_{k,m} \tag{3}$$

Due to the dense deployment of the SCs, some SCs can be turned off without affecting the QoS of the UEs. Thus, the total power consumption of SC *m* can be calculated as [13]:

$$P\_{m\_{\rm tot}} = \beta P\_{m\_{\rm on}} + (1 - \beta)P\_{m\_{\rm on}}\gamma\_m + \theta\_m P\_{m\_{\rm tx}}\gamma\_m \tag{4}$$

where *Pmon* and *Pmtx* are the baseline and the transmission power consumption, respectively, while *β* is the inactive level of the SC, such that *Pmo f f* = *βPmon* , and *θ<sup>m</sup>* is the portion of power consumption that is due to the feeder losses and power amplifier of SC *m* [51]. The total power consumption of the system (*Psys*) is the sum of the power consumption of all SCs. The PE of the system is given by [9]:

$$PE\_{sys} = \frac{\mathcal{C}\_{sys}}{P\_{sys}}\tag{5}$$

To improve the PE, the SCs on/off switching decisions using BPSO is proposed in this paper. The PSO is an iterative population-based search algorithm inspired by the hunting behavior of a flock of flying birds [52–54]. In PSO, every particle is considered a bird of the flock and represents a possible solution to the problem [30,55–57]. The search begins with an initial set of particles and attempts to find the best solution by searching around the solution space. The motion of the particle is based on its local best position, and the best-known position of all the other particles [9,52]. The fitness value of every particle is calculated using a fitness function that is optimized in every iteration [9,52].

The set of particles X = {*x*1, *x*2, . . . , *xNpar* } is defined, where *xNpar* represents one possible status for the SCs, while *Npar* is the swarm size. In BPSO algorithm, the population is randomly initialized as binary values. For every particle, the population binary value of 1 signifies the active SC, while 0 signifies the sleeping SC. The velocity of the particle *j* is initialized as [9,58]:

$$
\sigma\_j = \upsilon\_{\min} + \left(\upsilon\_{\max} - \upsilon\_{\min}\right) a\_1 \tag{6}
$$

where *vmin* and *vmax* denote the minimum and the maximum velocity of the particle, respectively, while *a*<sup>1</sup> is a random number uniformly distributed between 0 and 1 [9]. The

velocity and the position of the particles are updated in each iteration. The velocity of the particle *j* in iteration (*z* + 1) is updated as [9,26]:

$$\left(\boldsymbol{v}\_{j}\right)^{z+1} = \omega\_{z}\left.\boldsymbol{v}\_{j}\right|^{z} + c\_{1}\left.a\_{2}\left(\boldsymbol{P}\_{\text{best}\_{j}} - \boldsymbol{\chi}\_{j}\right^{z}\right) + c\_{2}\left.a\_{3}\left(\boldsymbol{G}\_{\text{best}} - \boldsymbol{\chi}\_{j}\right^{z}\right) \tag{7}$$

where *ω<sup>z</sup>* is the IW in the *z th* iteration, while *v z j* and *x z j* are the velocity and position of the particle *j* in the *z th* iteration, respectively. The best position of the particle *j* is denoted as *Pbest<sup>j</sup>* , while the global best position of all the particles is denoted as *Gbest*. Additionally, *c*<sup>1</sup> and *c*<sup>2</sup> denote the acceleration parameters [9]. Moreover, *a*<sup>2</sup> and *a*<sup>3</sup> are two random numbers uniformly distributed between 0 and 1. It is noted that the particle *x<sup>j</sup>* is a binary vector and the velocity *v<sup>j</sup>* is also a vector. The sigmoid function (*Sig vj*(*m*) *z* ) is given as [52,55]:

$$\operatorname{Sign}\left(v\_j(m)^z\right) = \frac{1}{1 + e^{-v\_j(m)^z}}\tag{8}$$

Thus, the on/off state of the SC *m* in particle *j* in the *z*-th iteration is calculated as [52,55]:

$$\text{tr}\_{\vec{\jmath}}(m)^z = \begin{cases} 1, & a\_4 < \text{Sign}\left(v\_{\vec{\jmath}}(m)^z\right) \\ 0, & \text{otherwise} \end{cases} \tag{9}$$

where *a*<sup>4</sup> is a random number uniformly distributed between 0 and 1.

Since the IW is the pivotal factor in the convergence of the PSO, it should be carefully adjusted. Thus, linearly increasing IW is utilized in this paper, where the IW linearly increases every iteration from *ωmin* to *ωmax* . The IW in the *z th* iteration is given as [33]:

$$
\omega\_z = (\omega\_{\text{max}} - \omega\_{\text{min}}) \left( \frac{z - 1}{Z\_{\text{max}} - 1} \right) + \omega\_{\text{min}} \tag{10}
$$

where *ωmax* and *ωmin* denote the maximum and minimum IW, respectively [33], and *Zmax* is the maximum number of iterations [59].

#### **4. Proposed Algorithms**

To alleviate the number of active SCs and enhance the PE of the system, a linearly increasing IW-BPSO algorithm for SC on/off switching is proposed in this paper. Moreover, a novel SFR technique using CTs is proposed for SC sub-band allocation. The BPSO algorithm is applied first while the linearly increasing IW enhances the convergence of the algorithm. Then, the sub-bands are allocated to the active SCs using the novel SFR technique. It is worth noting that the new SFR technique is applied after the SC switching to allocate the sub-bands to the active SCs, only aiming to reduce the number of operations needed in the sub-band allocation to the SCs

#### *4.1. SC on/off Switching Using Linearly Increasing IW-BPSO Algorithm*

In this paper, SC switching utilizing a linearly increasing IW-BPSO algorithm is proposed. The UE can associate with an SC *m* if *SINRk*,*<sup>m</sup>* exceeds a certain threshold (*SINRthr*). At the beginning of our proposed algorithm, each UE *k* calculates the SINR from all SCs to determine all SCs that it can possibly associate with, then sorts the received SINR from these SCs in a descending order. The UE associates with the SC with the highest received SINR. If the number of UEs in this SC is larger than the maximum number of UEs in the SC (*Umax*), the UE connects with the SC having the next highest SINR. This continues until every UE is associated with one SC.

To minimize the number of active SCs and to elevate the PE, it is required to switch off the excessive SCs, taking into consideration the QoS of the UEs, since the SINR of every UE exceeds *SINRthr*. Thus, a multi-objective optimization problem had to be solved; this problem can be written as:

$$\text{Objective1}: \min \left(\sum\_{m=1}^{M} \gamma\_m \right) \\ \text{Objective2}: \max \left(PE\_{\text{sys}} \right) \\ \tag{11}$$

subject to:

$$
\gamma\_{m\nu} \in \{0, 1\}, \forall m \in M \tag{11a}
$$

$$
\varphi\_{m,k} \in \{0, 1\}, \forall m \in M, \ k \in \mathcal{U}\_m \tag{11b}
$$

$$\sum\_{m=1}^{M} \, \_{m} \varphi\_{m,k} \, \gamma\_{m} = 1 \, \text{\textquotedblleft } m \in M \text{\textquotedblright} \text{\textquotedblleft } L \in \mathcal{U}\_{m} \tag{11c}$$

$$\mathcal{U}\_m \le \mathcal{U}\_{\text{max}} \; \forall \; m \in \mathcal{M} \tag{11d}$$

where constraint (11a) is the on/off state indicator of the SC. Constraint (11b) is the UE association indicator. Constraint (11c) states that the UE is associated with only one active SC. Constraint (11d) indicates that the number of UEs in an SC cannot exceed the maximum number of UEs in an SC (*Umax*).

Algorithm 1 summarizes the proposed BPSO-based on/off SC switching algorithm. In line 1, after the random binary initialization of the population of every particle, the velocity of the particles is initialized as Equation (6). Then, for each iteration till the maximum number of iterations (*Zmax*) is reached, the IW (*ωz*) is computed (line 3) and the position (*xj* ) and velocity (*v<sup>j</sup>* ) of every particle *j* are updated using Equations (7) and (9), respectively (lines 5–6). Next, the fitness value of every particle *j* (F (*x<sup>j</sup>* )) is computed using Equation (11) (line 7). If *F xj* < *F Pbest<sup>j</sup>* , then *Pbest<sup>j</sup>* = *x<sup>j</sup>* , and if *F xj* < *F*(*Gbest*) , then *Gbest* = *x<sup>j</sup>* (lines 8–13).

**Algorithm 1:** Proposed linearly increasing IW-BPSO-based on/off SC switching.

**Inputs:** Locations of UEs, locations of SCs, swarm size (*Npar*), maximum number of iterations (*Zmax* ) **Output:** SC on/off indicator

1: **Initialize** the position (*x<sup>j</sup>* ) **randomly** and velocity (*v<sup>j</sup>* ) of every particle *j* as Equation (6).

2: **For** *z* = 1 to *Zmax* 3: **Calculate** *ωz* using Equation (10) 4: **For** each particle *j* 5: **Update** *v<sup>j</sup>* using Equation (7) 6: **Update** *x<sup>j</sup>* using Equation (9) 7: **Calculate** new fitness value *F* (*x<sup>j</sup>* ) as Equation (5) 8: **if** *F xj* < *F Pbest<sup>j</sup>* 

9: *Pbest<sup>j</sup>* ← *x<sup>j</sup>* 10: **end if** 11: **if** *F xj* < *F*(*Gbest*) 12: *Gbest* ← *x<sup>j</sup>* 13: **end if** 14: **end For** 15: **end For**

#### *4.2. SC Sub-Band Allotment Using Classification Trees (CTs)*

After the on/off switching decisions for all SCs are taken, the second phase is initialized, which is the sub-band allocation for the active SCs based on the SFR, which is illustrated in Figure 2, over three cells without loss of generality. In the shown example, if we have three hexagonal-shaped cells, each split into center and edge regions, the frequency band is divided to three (*Nsub* = 3) sub-bands: X, Y, and Z. The edge region of cells 1, 2, and 3 are allocated sub-bands X, Y, and Z, respectively. Consequently, cell 1 is allocated sub-bands Y and Z in the center region. Similarly, cell 2 is allocated X and Z, and cell 3 is allocated X and Y.

bands Y and Z in the center region. Similarly, cell 2 is allocated X and Z, and cell 3 is

**Figure 2.** SFR example for hexagonal shaped cells with = 3. **Figure 2.** SFR example for hexagonal shaped cells with *Nsub* = 3.

allocated X and Y.

In this technique, every SC is divided into center and edge regions. One of the sub-bands can be used in the edge region of every SC, on condition that it is not used by the edge regions of the adjoin SCs. The center region of every SC can use the remaining sub-bands with reduced transmission power. This alleviates not only the interference to the adjoin SCs but also the power consumption of the whole network. To determine the vertices of the center region, the distance between the center of the SC and every vertex of the SC is computed. Then the distance between the center of the SC and the nearest SC vertex (the smallest distance) is determined and is regarded as the SC radius. The radius of the center zone is chosen as 50% of the SC radius as it maximizes the throughput of the system [13]. A real example is demonstrated in Figure 3a, displaying an SC (the purple In this technique, every SC is divided into center and edge regions. One of the *Nsub* sub-bands can be used in the edge region of every SC, on condition that it is not used by the edge regions of the adjoin SCs. The center region of every SC can use the remaining sub-bands with reduced transmission power. This alleviates not only the interference to the adjoin SCs but also the power consumption of the whole network. To determine the vertices of the center region, the distance between the center of the SC and every vertex of the SC is computed. Then the distance between the center of the SC and the nearest SC vertex (the smallest distance) is determined and is regarded as the SC radius. The radius of the center zone is chosen as 50% of the SC radius as it maximizes the throughput of the system [13]. A real example is demonstrated in Figure 3a, displaying an SC (the purple SC) and its adjoining SCs. Figure 3b displays the seven sub-bands. Every SC uses one of the seven sub-bands in its edge region. While the center region of this SC (the grey region) can use the remaining six sub-bands.

SC) and its adjoining SCs. Figure 3b displays the seven sub-bands. Every SC uses one of the seven sub-bands in its edge region. While the center region of this SC (the grey region) can use the remaining six sub-bands. First, the SC senses the signals of the edge region of the adjoin SCs upon switching to obviate utilizing them. A binary indicator, *τ<sup>m</sup>* = 1, indicates the presence of unused sub-bands by the edge region of the adjoin SCs; otherwise, *τ<sup>m</sup>* = 0. If there are unused sub-bands by the edge region of the adjoined SCs (*τ<sup>m</sup>* = 1), then for every unused sub-band *f*, the distance between the SC and the closest nonadjacent SC using sub-band *f* in its edge region (*DNon*−*ad f*) is measured, and this is repeated for the rest of the remaining unused sub-bands. If all the sub-bands are used (*τ<sup>m</sup>* = 0), then for every sub-band *q,* the distance between the SC and the adjoin SC using sub-band *q* in its edge region (*DAdq*) is measured. Then, the edge region of the SC is allocated the sub-band used in the farthest adjoined SC.

(**a**) (**b**)

**Figure 3.** (**a**) A real example demonstrating a SC (the purple SC) with the center region (the grey region); (**b**) the seven used sub-bands. **Figure 3.** (**a**) A real example demonstrating a SC (the purple SC) with the center region (the grey region); (**b**) the seven used sub-bands.

First, the SC senses the signals of the edge region of the adjoin SCs upon switching to obviate utilizing them. A binary indicator, = 1, indicates the presence of unused sub-bands by the edge region of the adjoin SCs; otherwise, = 0. If there are unused sub-bands by the edge region of the adjoined SCs ( = 1), then for every unused subband *f*, the distance between the SC and the closest nonadjacent SC using sub-band *f* in its edge region (ேିௗ) is measured, and this is repeated for the rest of the remaining unused sub-bands. If all the sub-bands are used ( = 0), then for every sub-band *q,* the distance between the SC and the adjoin SC using sub-band *q* in its edge region (ௗ) is measured. Then, the edge region of the SC is allocated the sub-band used in the farthest According to *<sup>τ</sup><sup>m</sup>* and *<sup>D</sup>Non*−*ad f* or *<sup>D</sup>Adq*, the CT takes the decision to assign which sub-band to the edge region of the SC. A CT example of 3 sub-bands is shown in Figure 4. In case of the presence of unused sub-bands in the edge region of the adjoin SCs (*τ<sup>m</sup>* = 1), then if the distance between the SC and the closest SC using the sub-band *1* in its edge region (*DNon*−*ad f* = *<sup>D</sup>Non*−*ad*1) is larger than *<sup>D</sup>Non*−*ad*2; *<sup>D</sup>Non*−*ad*<sup>1</sup> and *<sup>D</sup>Non*−*ad*<sup>3</sup> are checked. If *DNon*−*ad*<sup>1</sup> is larger than *DNon*−*ad*3, then sub-band *<sup>1</sup>* is selected. In case all the sub-bands are utilized (*τ<sup>m</sup>* = 0), then if the distance between the SC and the adjacent SC using sub-band *2* (*DAdq* = *DAd*2) is larger than *DAd*1, then *DAd*<sup>2</sup> and *DAd*<sup>3</sup> are checked. If *DAd*<sup>3</sup> is larger than *DAd*2, then sub-band *3* is selected. Afterward, the center region can utilize the remaining sub-bands with lower transmission power.

adjoined SC. According to and ேିௗ or ௗ, the CT takes the decision to assign which sub-band to the edge region of the SC. A CT example of 3 sub-bands is shown in Figure 4. In case of the presence of unused sub-bands in the edge region of the adjoin SCs ( = 1), then if the distance between the SC and the closest SC using the sub-band *1* in its edge region ( ேିௗ = ேିௗଵ) is larger than ேିௗଶ ; ேିௗଵ and ேିௗଷ are checked. If ேିௗଵ is larger than ேିௗଷ, then sub-band *1* is selected. In case all the sub-bands are utilized ( = 0), then if the distance between the SC and the adjacent SC using sub-band *2* (ௗ = ௗଶ) is larger than ௗଵ, then ௗଶ and ௗଷ are checked. If ௗଷ is larger than ௗଶ, then sub-band *3* is selected. Afterward, the center region can utilize the remaining sub-bands with lower transmission power. Unlike the computational complexity of the conventional SFR algorithm *O*(*S* 2 *. Nsub* 2 ), which depends on the network size (*S*) and the number of sub-bands (*Nsub*) [60], the computational complexity of the proposed CTs algorithm is much lower. Since the computational complexity of the DTs is *O(1)* [61], as no multiplication process is done and only a sequence of branching operations are performed, then computational complexity of the CTs is also *O(1)* for every SC. Thus, the computational complexity of the system is *S* x *O(1),* since only a sequence of branching operations are performed on moving along the CT according to the binary indicator *(τm*) and the distance between the SC and the closest nonadjacent SC (*DNon*−*ad f*)/farthest adjoining SC (*DAd<sup>q</sup>* ) using the sub-band in its edge region. Hence, applying SFR after the SC switching greatly reduces the computational complexity as SFR will be applied to the active SCs only (smaller *S*), while the complexity of the BPSO algorithm *O(Zmax . P)* depends on the maximum number of iterations and the population size (*P*) [62,63]. On the other hand, the average computational time of the proposed algorithm is less than one minute (54.857 s). In modern HetNets, the traffic load of the network is monitored for longer time periods in the order of 5–15 min [6,64–66], while the on-off switching decisions of the SCs are usually taken every 15–60 min [64,66]. Thus, the proposed algorithm is suitable for practical implementation in real time.

**Figure 3.** (**a**) A real example demonstrating a SC (the purple SC) with the center region (the grey

used sub-bands. If all the sub-bands are used ( = 0), then for every sub-band *q,* the distance between the SC and the adjoin SC using sub-band *q* in its edge region (

measured. Then, the edge region of the SC is allocated the sub-band used in the farthest

sub-band to the edge region of the SC. A CT example of 3 sub-bands is shown in Figure 4. In case of the presence of unused sub-bands in the edge region of the adjoin SCs ( = 1), then if the distance between the SC and the closest SC using the sub-band *1* in its edge

) is larger than −2

) is larger than 1

sub-bands are utilized ( = 0), then if the distance between the SC and the adjacent SC

is larger than −3

utilize the remaining sub-bands with lower transmission power.

or

First, the SC senses the signals of the edge region of the adjoin SCs upon switching to obviate utilizing them. A binary indicator, = 1, indicates the presence of unused sub-bands by the edge region of the adjoin SCs; otherwise, = 0. If there are unused sub-bands by the edge region of the adjoined SCs ( = 1), then for every unused subband *f*, the distance between the SC and the closest nonadjacent SC using sub-band *f* in its

) is measured, and this is repeated for the rest of the remaining un-

, the CT takes the decision to assign which

, then sub-band *1* is selected. In case all the

and 3

and −3

are checked. If

; −1

, then 2

, then sub-band *3* is selected. Afterward, the center region can

) is

are

region); (**b**) the seven used sub-bands.

According to and −

region ( − = −1

using sub-band *2* ( = 2

is larger than 2

checked. If −1

edge region (−

adjoined SC.

3

**Figure 4.** Classification tree (CT) example for three sub-bands. **Figure 4.** Classification tree (CT) example for three sub-bands.

#### Unlike the computational complexity of the conventional SFR algorithm *O*(*S². ²*), **5. Numerical Results**

which depends on the network size (*S*) and the number of sub-bands () [60], the computational complexity of the proposed CTs algorithm is much lower. Since the computational complexity of the DTs is *O(1)* [61], as no multiplication process is done and only a The simulation parameters are presented in Table 1. Various Voronoi SCs are allocated in the coverage area of the MCs. The UEs are randomly deployed within the entire network. The simulations are performed utilizing MATLAB R2018. The swarm size is chosen in the range of 20–50 particles [9,67–71]. It is shown in [71] that 25 is the optimum size over 16 different sizes on 3 different cases. The IW is linearly increased from 0.4 to 0.9, which is the range recommended in [33,72,73]. The velocity of the particles is chosen in the range of [−0.6, 0.6] [9,74,75]. Adjusting these parameters is essential for the convergence of the BPSO algorithm.

The performance of several algorithms is assessed, they are as follow:


The number of active SCs for various number of UEs is demonstrated in Figure 5. The number of active SCs increases with increasing the number of UEs, as more SCs are being activated to guarantee the minimum required SINR of the UEs allowing the UEs to associate with the SC having the better SINR, enhancing their data rates and improving the system performance. The "*ω* = 0.4" has better convergence than the "*ω* = 0.9" for the 200, 500, and 900 UEs, respectively, since reducing the IW facilitates in exploring the search space (global search), while raising the IW aids in exploiting (local search) the search space. The "linearly increasing *ω*" shows better convergence than the fixed "*ω* = 0.9" and "*ω* = 0.4", as in the "linearly increasing *ω*" case, the IW is lower at the beginning, allowing better exploration of the search space, and then it linearly increases, enhancing the local search. The fixed IW can get trapped during the search and is unable to find the global minimum number of active SCs. Therefore, the linearly increasing IW is chosen to be utilized in the proposed algorithm for the rest of the simulations in this paper.


be utilized in the proposed algorithm for the rest of the simulations in this paper.

**Figure 5.** Number of active SCs for various numbers of UEs for different values of IW. **Figure 5.** Number of active SCs for various numbers of UEs for different values of IW.

UEs, since more SCs are switched on to sustain the minimum SINR of the UEs.

The number of active SCs for various numbers of UEs is presented in Figure 6. The "always on" algorithm has the highest number of active SCs, since all the SCs are kept active. The "random 10%" algorithm has a lower number of active SCs than the "always on" algorithm, since in the "random 10%" algorithm, 10% of the SCs are randomly turned off. However, the number of active SCs is constant with regard to the number of UEs. The proposed and "BPSO only" algorithms have the least number of active SCs, as the SCs are turned off utilizing the BPSO algorithm. In both algorithms (BPSO only and proposed), BPSO is used for SC switching. That is why the number of active SCs of the two algorithms is very close. Moreover, the number of active SCs increases on increasing the number of The number of active SCs for various numbers of UEs is presented in Figure 6. The "always on" algorithm has the highest number of active SCs, since all the SCs are kept active. The "random 10%" algorithm has a lower number of active SCs than the "always on" algorithm, since in the "random 10%" algorithm, 10% of the SCs are randomly turned off. However, the number of active SCs is constant with regard to the number of UEs. The proposed and "BPSO only" algorithms have the least number of active SCs, as the SCs are turned off utilizing the BPSO algorithm. In both algorithms (BPSO only and proposed), BPSO is used for SC switching. That is why the number of active SCs of the two algorithms

The total system throughput for various numbers of UEs is shown in Figure 7. The proposed algorithm has the largest system throughput, since the proposed algorithm alleviates the interference levels, as it utilizes the SFR and switches the SCs on/off using the

*Sensors* **2022**, *22*, x FOR PEER REVIEW 12 of 21

is very close. Moreover, the number of active SCs increases on increasing the number of UEs, since more SCs are switched on to sustain the minimum SINR of the UEs. cause of the random selection of the switched-off SCs, which is not the optimum one. The "always on" algorithm has the minimum system throughput because it has the largest

interference levels as all the SCs are continuously active and it does not utilize the SFR.

proposed algorithm because of its larger interference levels for not using SFR. The "random 10%" algorithm has lower system throughput than the "BPSO only" algorithm be-

**Figure 6.** Number of active SCs for various numbers of UEs in case of linearly increasing IW. **Figure 6.** Number of active SCs for various numbers of UEs in case of linearly increasing IW.

The total system throughput for various numbers of UEs is shown in Figure 7. The proposed algorithm has the largest system throughput, since the proposed algorithm alleviates the interference levels, as it utilizes the SFR and switches the SCs on/off using the BPSO algorithm. The "BPSO only" algorithm has lower system throughput than the proposed algorithm because of its larger interference levels for not using SFR. The "random 10%" algorithm has lower system throughput than the "BPSO only" algorithm because of the random selection of the switched-off SCs, which is not the optimum one. The "always on" algorithm has the minimum system throughput because it has the largest interference levels as all the SCs are continuously active and it does not utilize the SFR.

The total power consumption of the system for various numbers of UEs is presented in Figure 8. The "always on" algorithm has the largest power consumption since no switchingoff techniques are used. The "random 10%" algorithm has less power consumption than the "always on" algorithm because of the power savings from the switched-off SCs and higher power consumption than the "BPSO only" algorithm. The proposed algorithm consumes the least power because of the turning-off of the SCs utilizing BPSO, then applying the SFR minimizing the power consumption furthermore, because of the reduced transmission power for the UEs in the center region. In the proposed algorithm, the power consumption increases with the increment of the number of UEs, as more SCs are kept active to ensure the QoS of the UEs.

*Sensors* **2022**, *22*, x FOR PEER REVIEW 13 of 21

**Figure 7.** Total system throughput for various numbers of UEs. **Figure 7.** Total system throughput for various numbers of UEs.

*Sensors* **2022**, *22*, x FOR PEER REVIEW 14 of 21

The total power consumption of the system for various numbers of UEs is presented

The PE of the system for various numbers of UEs is demonstrated in Figure 9. The proposed algorithm has the highest PE because of the reduced power consumption and the elevated system throughput. The "BPSO only" algorithm has lower PE because of its alleviated throughput and higher system power consumption compared to the proposed algorithm. The "random 10%" has lower PE than the "BPSO only" algorithm because of its alleviated system throughput caused by the random choice of the switched-off SCs. The

The outage probability for various SINR thresholds in the case of 900 UEs is depicted in Figure 10. The outage probability is defined as the percentage of UEs that are unable to attain a certain SINR threshold. The proposed algorithm has the minimum outage probability because of the enhanced SINR of the users, which resulted from minimizing the interference. The "BPSO only" algorithm has larger outage probability than the proposed algorithm, since it has larger interference levels. Both algorithms consider the QoS of the

10%" algorithm has larger outage probability than the "BPSO only" because of the deteriorated system performance. The "always on" algorithm has the largest outage probability because of the huge interference levels, as the SFR principle is not utilized and no on/off

switching techniques are used, diminishing the total performance of the system.

in both algorithms. The "random

**Figure 8.** Total system power consumption for various numbers of UEs. **Figure 8.** Total system power consumption for various numbers of UEs.

users, as all the users have SINR larger than ℎ

the least data rates because all the SCs are continuously active.

The PE of the system for various numbers of UEs is demonstrated in Figure 9. The proposed algorithm has the highest PE because of the reduced power consumption and the elevated system throughput. The "BPSO only" algorithm has lower PE because of its alleviated throughput and higher system power consumption compared to the proposed algorithm. The "random 10%" has lower PE than the "BPSO only" algorithm because of its alleviated system throughput caused by the random choice of the switched-off SCs. The "always on" algorithm has the least PE, since it has the largest power consumption and the least data rates because all the SCs are continuously active.

The outage probability for various SINR thresholds in the case of 900 UEs is depicted in Figure 10. The outage probability is defined as the percentage of UEs that are unable to attain a certain SINR threshold. The proposed algorithm has the minimum outage probability because of the enhanced SINR of the users, which resulted from minimizing the interference. The "BPSO only" algorithm has larger outage probability than the proposed algorithm, since it has larger interference levels. Both algorithms consider the QoS of the users, as all the users have SINR larger than *SINRthr* in both algorithms. The "random 10%" algorithm has larger outage probability than the "BPSO only" because of the deteriorated system performance. The "always on" algorithm has the largest outage probability because of the huge interference levels, as the SFR principle is not utilized and no on/off switching techniques are used, diminishing the total performance of the system. *Sensors* **2022**, *22*, x FOR PEER REVIEW 15 of 21

**Figure 9.** Power efficiency for various numbers of UEs. **Figure 9.** Power efficiency for various numbers of UEs.

**Figure 10.** Outage probability for various SINR thresholds in the case of 900 UEs.

**6. Conclusions**

**Figure 9.** Power efficiency for various numbers of UEs.

**Figure 10.** Outage probability for various SINR thresholds in the case of 900 UEs. **Figure 10.** Outage probability for various SINR thresholds in the case of 900 UEs.

#### **6. Conclusions 6. Conclusions**

Minimizing the power consumption of the SCs and elevating the PE of the network are huge challenges facing the 5G HetNets. In this paper, to tackle these challenges, novel algorithms are proposed based on linear increasing IW-BPSO and SFR. The BPSO algorithm is used for SC on/off switching reducing the power consumption of the system without deteriorating the QoS of the UEs. Moreover, the linearly increasing IW is exploited to enhance the convergence of the BPSO algorithm to find the minimum number of active SCs. Furthermore, the CT-based SFR is proposed, where the SCs are divided into center and edge regions and different sub-bands are allocated to the edge regions of the adjoining SCs, minimizing the interference among the SCs. The results demonstrate that the proposed algorithms surpass the other conventional algorithms with regard to the total system power consumption, the total system throughput, the PE, and the outage probability. Additional work can be done in the future to address the coverage hole problem in Voronoi cells and to enhance the accuracy of the PSO algorithm using dynamic inertia weight.

**Author Contributions:** Conceptualization, M.O., B.A. and S.E.R.; methodology, M.O., B.A. and S.E.R.; software M.O.; writing—original draft preparation, M.O.; writing—review and editing, B.A. and S.E.R.; All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**


#### **References**

