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

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.

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. 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:

• User scheduling: Due to that only a few users can be bound into a DVB-S2X frame, and there are a large number of active users in each multicast group, it is necessary to design the user scheduling algorithm. Meanwhile, it should be noted that the interference between scheduled users depends on the beamforming design, which in turn depends on the scheduled users in other beams. Therefore, user scheduling and beamforming design are coupled, and the joint design scheme of user scheduling and beamforming needs to be considered.
• Channel errors: The LEO satellite has high orbital speed, which will produce a large Doppler shift and result in the channel phase deviation [8]. Meanwhile, the factors such as distortion of high-frequency devices, expiration of the CSI and large propagation delay also can cause the channel phase disturbance. Therefore, it is difficult to obtain accurate channel state information (CSI) at the LEO satellite transmitter. Due to the existence of CSI errors, the designed beamforming vector does not match the actual CSI, resulting in the reduction in the receiving gain and signal to interference plus noise ratio (SINR) of the user terminal. Then, the QoS will not be guaranteed. Thus, it is of practical significance to study the robust user scheduling and beamforming design.
• EE optimization: Due to the limited energy load of LEO satellites, to prolong the service life of the LEO satellite and improve the stability of the LEO satellite communications system, under the consideration of green communications and economic benefits, we need to pay attention to the EE optimization [9].
• Beamforming scheme: In the beamforming architectures of the massive MIMO technology, although the digital beamforming design can significantly improve the SE, it would bring high hardware complexity and high power consumption. Although the hardware overhead of hybrid beamforming architecture based on full connection is slightly higher than that of the partial connection architecture, it can balance the hardware complexity and system performance, and has higher cost performance. Therefore, in this paper, we selected the hybrid beamforming technology based on the full connection structure [10].

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:

• We establish the downlink transmission system model and channel model of the massive MIMO LEO satellite multigroup multicast communication system and analyze the CSI errors.
• Based on the CSI, we adopt a low complexity hierarchical clustering algorithm based on the $Ward$ connection method to group users, which can lay a foundation for the joint user scheduling and beamforming design.
• We establish the joint user scheduling and hybrid beamforming design problem model based on EE maximization, and binary variables are defined to represent whether the user is scheduled or not. Then, we transform the optimization problem into a Boolean fractional programming (BFP) problem, which is also a quadratic constraint quadratic programming (QCQP) form problem.
• For the BFP problem in QCQP form, we invoke the quadratic transformation algorithm to handle the fractional programming form problem in the objective function. Meanwhile, the SDP algorithm is invoked to convert the objective function in QCQP 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 transformed into a difference of convex (DC) programming problem.
• 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 SDP algorithm, a penalty iterative algorithm is adopted.
• For the solution of the digital beamforming matrix and the analog beamforming matrix in the hybrid beamformer, the MM-AltOpt algorithm is proposed.

3. System Model and Problem Formulation

3.1. System Model

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_x\times N_y$ antennas and $L\left(L\le N\right)$ RF links and covers $L$ multicast groups and $K$ active users, where $K\ge L,K\ge N$. Let the multicast group covered by the $l\mathrm{th}$ beam be $U_l$ and the number of users in this multicast group be $\left|U_l\right|$, assuming that the number of users that can be accommodated in each DVB-S2X frame is $U_s$. It should be noted that each user terminal belongs to only one multicast group, i.e., $U_i\cap U_j=\mathsf{\varnothing },\forall i,j\in \left\{1,\dots ,L\right\},i\ne 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.

The received signal $y_{k,l}$ of the $k\mathrm{th}$ user in the $l\mathrm{th}$ multicast group can be expressed as:

$$y_{k,l}={\mathit{h}}_{k,l}^H{\mathit{F}}_{RF}{\mathit{F}}_{BB}[:,l]s_l+{\displaystyle \sum _{j=1,j\ne l}^L{\mathit{h}}_{k,j}^H{\mathit{F}}_{RF}{\mathit{F}}_{BB}[:,j]}{s}_l+n_{k,l},k\in \left\{1,\dots ,\left|U_l\right|\right\},l\in \left\{1,\dots ,L\right\},$$

(1)

where the first term in (1) represents the expected received signal of the $k\mathrm{th}$ user in the $l\mathrm{th}$ multicast group, the second term represents the interference of other multicast groups and the third term represents the additive Gaussian white noise; ${\mathit{h}}_{k,l}\in {ℂ}^{N\times 1}$ represents the channel vector of the $k\mathrm{th}$ user in the $l\mathrm{th}$ multicast group, ${\mathit{F}}_{BB}\in {ℂ}^{L\times L}$ represents the digital beamforming matrix, ${\mathit{F}}_{BB}[:,l]\in {ℂ}^{L\times 1}$ represents the digital beamforming vector of the $l\mathrm{th}$ multicast group, and ${\mathit{F}}_{RF}\in {ℂ}^{N\times L}$ represents the analog beamforming matrix, where each element of ${\mathit{F}}_{RF}$ should meet the unit modulus element [19], i.e., $\left|{\left({\mathit{F}}_{RF}\right)}_{i,j}\right|=1$. In addition, $s_l$ represents the signal of the multicast group $U_l$, which meets the unit power constraint, i.e., ${\rm E}\left\{{\left|s_l\right|}^2\right\}=1$, $n_l\sim CN(0,{\sigma }^2)$ represents the additive Gaussian white noise, which is related to the Boltzmann constant $\kappa$, system bandwidth $B$ and the noise temperature $T$.

For the convenience of analysis, we set $\mathit{F}\in {ℂ}^{N\times L}={\mathit{F}}_{RF}{\mathit{F}}_{BB}=[{\mathit{f}}_1,{\mathit{f}}_2,\dots ,{\mathit{f}}_L]$ as the hybrid beamforming matrix, and ${\mathit{f}}_l\in {ℂ}^{N\times 1}={\mathit{F}}_{RF}{\mathit{F}}_{BB}[:,l]$ is the hybrid beamforming vector of the $l\mathrm{th}$ multicast group. Therefore, (1) can be rewritten as:

$$y_{k,l}={\mathit{h}}_{k,l}^H{\mathit{f}}_l{s}_l+{\displaystyle \sum _{j=1,j\ne l}^L{\mathit{h}}_{k,j}^H{\mathit{f}}_j}{s}_l+n_{k,l},k\in \left\{1,\dots ,\left|U_l\right|\right\},l\in \left\{1,\dots ,L\right\}$$

(2)

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\mathrm{th}$ user in the $l\mathrm{th}$ multicast group at instant $t$ and frequency $f$ can be modeled as follows [20]:

$${\mathit{h}}_{k,l}(t,f)={\displaystyle \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 {\mathit{V}}_{k,l,p},$$

(3)

where $f$ denotes the carrier frequency, $a_{k,l,p}$,$f_{d(k,l,p)}$,${\tau }_{k,l,p}$ are the complex channel gain, Doppler shift and propagation delay, respectively, $P_{k,l}$ denotes the number of propagation paths and ${\mathit{V}}_{k,l,p}\in {ℂ}^{N\times 1}$ is the UPA array response vector, which can be given by

$${\mathit{V}}_{k,l,p}=\mathit{V}({\phi }_{k,l,p}^x,{\phi }_{k,l,p}^y)={\mathit{v}}_{N_x}\left(\sin {\phi }_{k,l,p}^y\cos {\phi }_{k,l,p}^x\right)\otimes {\mathit{v}}_{N_y}\left(\cos {\phi }_{k,l,p}^y\right),$$

(4)

$${\mathit{v}}_{N_x}(\sin {\phi }_{k,l,p}^y\cos {\phi }_{k,l,p}^x)=\frac{1}{\sqrt{N_x}}\left(1,e^{-j\frac{2\pi d}{\lambda }\sin {\phi }_{k,l,p}^y\cos {\phi }_{k,l,p}^x},\dots ,e^{-j\frac{2\pi d}{\lambda }\left(N_x-1\right)\sin {\phi }_{k,l,p}^y\cos {\phi }_{k,l,p}^x}\right),$$

(5)

$${\mathit{v}}_{N_y}(\cos {\phi }_{k,l,p}^y)=\frac{1}{\sqrt{N_y}}\left(1,e^{-j\frac{2\pi d}{\lambda }\cos {\phi }_{k,l,p}^y},\dots ,e^{-j\frac{2\pi d}{\lambda }\left(N_y-1\right)\cos {\phi }_{k,l,p}^y}\right),$$

(6)

where ${\phi }_{k,l,p}^x$, ${\phi }_{k,l,p}^y$ represent azimuth angle and pitch angle associated with the propagation path $p$ of the $k\mathrm{th}$ user in the $l\mathrm{th}$ multicast group, respectively, $\lambda$ denotes the wavelength and $d$ represents the spacing of antenna elements, the value of which is usually $\lambda /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:

$${\mathit{h}}_{k,l}={\overline{\mathit{h}}}_{k,l}+{\tilde{\mathit{h}}}_{k,l},$$

(7)

where ${\overline{\mathit{h}}}_{k,l}=\sqrt{\frac{{\kappa }_{k,l}{\gamma }_{k,l}}{{\kappa }_{k,l}+1}}\times {\mathit{V}}_{k,l}$ represents the LoS component, ${\tilde{\mathit{h}}}_{k,l}=\sqrt{\frac{{\gamma }_{k,l}}{{\kappa }_{k,l}+1}}\times {\mathit{V}}_{k,l,c}\times {\mathit{V}}_{k,l}^H$ represents the multipath component, ${\kappa }_{k,l}$ denotes the Rician factor, ${\mathit{V}}_{k,l,c}\in {ℂ}^{N_{ut}\times 1}\sim ℂ\mathbb{N}\left(0,\sum \right)$ represents Rician component, $Tr(\sum )=1$ and ${\gamma }_{k,l}$ represents the average channel power, which mainly includes the transmit antenna gain $G_{leo}$, receiver antenna gain $G_{ut}$ and link power loss. The link power loss is mainly caused by the free space path loss $L{P}_{fs}$ and the atmospheric absorption loss $L{P}_{at}$. Therefore, the average channel power ${\gamma }_{k,l}$ can be written as

$${\gamma }_{k,l}=G_{leo}[\mathrm{d}\mathrm{B}]+G_{ut}[\mathrm{d}\mathrm{B}]-L{P}_{at}[\mathrm{d}\mathrm{B}]-L{P}_{fs}[\mathrm{d}\mathrm{B}],$$

(8)

where $L{P}_{fs}$ can be given by $L{P}_{fs}=20\left({\log }_{10}\left(D_{k,l}\right)\right)+{\log }_{10}\left(f\right)+{\log }_{10}\left(4\pi /c\right)$, $c$ is the speed of light, $D_{k,l}$ represents the transmission distance, $L{P}_{at}$ is related to the carrier frequency, temperature $T\left(h\right)$, pressure $P\left(h\right)$ and humidity $\rho \left(h\right)$, which can be given by $L{P}_{at}={\displaystyle {\int }_{h_{ut}}^{h_{at}}L{P}_{at}\left(f,T\left(h\right),P\left(h\right),\rho \left(h\right)\right)dh}$, $L{P}_{at}\left(f,T\left(h\right),P\left(h\right),\rho \left(h\right)\right)$ is the loss per meter, $h_{ut}$ is the user’s height and $h_{at}$ is the atmosphere thickness. The specific calculation method of $L{P}_{at}$ can be found in the literature [21].

For the convenience of analysis, it is assumed that the parameters in the channel vector ${\mathit{h}}_{k,l}$ are constant within coherence time and change over time in a certain ergodic process. In (3), the Doppler shift $f_{d(k,l,p)}$ and the propagation delay ${\tau }_{k,l,p}$ 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_{d(k,l,p)}^{leo}$ generated by the LEO satellite motion and the Doppler shift $f_{d(k,l,p)}^{ut}$ generated by the users’ motion [22]. Since the transmission between the LEO satellite and user terminals is mainly under LOS, $f_{d(k,l,p)}^{leo}$ of different transmission paths can be considered to be the same, and we omit the path index of $f_{d(k,l,p)}^{leo}$, i.e., ${\left\{f_{d(k,l,p)}^{leo}\right\}}_1^{P_k}=f_{d(k,l)}^{leo}$; $f_{d(k,l)}^{leo}$ 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)}^{leo}=-\frac{f}{c}\times \frac{{\omega }^{leo}{r}_{e}r\sin \left({\varphi }_t-{\varphi }_{t_0}\right)\mu \left({\theta }_{\max }\right)}{\sqrt{r_e^2+r^2-2r_{e}r\cos \left({\varphi }_t-{\varphi }_{t_0}\right)\mu \left({\theta }_{\max }\right)}},$$

(9)

where $\mu \left({\theta }_{\max }\right)=\cos \left[{\cos }^{-1}\left(\frac{r_e}{r}\cos {\theta }_{\max }\right)-{\theta }_{\max }\right]$, $r_e$ denotes the earth radius, $r$ represents the distance between the LEO satellite track point and the earth center, $\left({\varphi }_t-{\varphi }_{t_0}\right)$ represents the angular distance of the earth’s surface along the LEO satellite trajectory from instant $t$ to instant $t_0$, $\omega$ represents the angular velocity of the LEO satellite and $c$ represents the speed of light.

The value of $f_{d(k,l,p)}^{ut}$ in different transmission paths is different, which is mainly caused by the movement of user terminals and surrounding scatterers. The power spectrum of $f_{d(k,l,p)}^{ut}$ 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-{(\frac{f_{d(k,l,p)}^{ut}}{f_{d(k,l,p),\max }^{ut}})}^2}}$$

(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_d^e$, the Doppler shift compensation of size $f_d^{cps}=f_d^e$ 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_d^{rsd}$, 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^e}^2=CRLB(f)=\frac{1}{\mathrm{S}\mathrm{N}\mathrm{R}}\frac{3}{2{\pi }^2{T}^{2}N(N^2-1)},$$

(11)

where $N$ is the pilot length, $T$ is the sampling time and $\mathrm{S}\mathrm{N}\mathrm{R}$ is the signal-to-noise ratio. According to the properties of variance, after Doppler shift compensation of size $f_d^{cps}$, the variance of residual Doppler shift $f_d^{rsd}$ is equal to the variance of $f_d^e$, i.e.,

$${\sigma }_{f_d^{rsd}}^2={\sigma }_{f_d^e}^2=\frac{1}{\mathrm{S}\mathrm{N}\mathrm{R}}\frac{3}{2{\pi }^2{T}^{2}N(N^2-1)},$$

(12)

The influence of residual Doppler shift on channel phase errors can be expressed as ${\varphi }_{f_d^{rsd}}=2\pi f_d^{rsd}\Delta T$, where $\Delta T$ represents the sum of the downlink propagation delay and the duration of the DVB-S2X frame. Therefore, the variance of ${\varphi }_{f_d^{rsd}}$ can be expressed as ${\sigma }_{{\varphi }_{f_d^{rsd}}}^2=4{\pi }^2{\sigma }_{f_d^{rsd}}^2\Delta T^2=\frac{1}{SNR}\frac{6\Delta T^2}{T^{2}N(N^2-1)}$, and ${\varphi }_{f_d^{rsd}}$ follows a real-valued Gaussian distribution with mean zero and variance ${\sigma }_{{\varphi }_{f_d^{rsd}}}^2$, i.e., ${\varphi }_{f_d^{rsd}}\sim N\left(0,{\sigma }_{{\varphi }_{f_d^{rsd}}}^2\right)$. The channel error vector caused by ${\varphi }_{f_d^{rsd}}$ can be expressed as ${\mathit{v}}_{{\varphi }_{f_d^{rsd}}}={[e^{j{\varphi }_{f_d^{rsd},1}},e^{j{\varphi }_{f_d^{rsd},2}},\dots ,e^{j{\varphi }_{f_d^{rsd},N}}]}^T$.

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 ${\tau }^{cps}=\beta {\tau }_{\min }+(1-\beta ){\tau }_{\max }, (0\le \beta \le 1)$ is usually implemented at the receivers, and ${\tau }_{k,l}^{\min }=\min {\left\{{\tau }_{k,l,p}\right\}}_1^{P_{k,l}}$ and ${\tau }_{k,l}^{\max }=\max {\left\{{\tau }_{k,l,p}\right\}}_1^{P_{k,l}}$ represent the minimum propagation delay and the maximum propagation delay of the $k\mathrm{th}$ user in the $l\mathrm{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 ${\varphi }_{\tau }$, which follows a real-valued Gaussian distribution with mean zero and variance ${\sigma }_{{\varphi }_{\tau }}^2$, i.e., ${\varphi }_{\tau }\sim N\left(0,{\sigma }_{{\varphi }_{\tau }}^2\right)$. The channel error vector caused by ${\varphi }_{\tau }$ can be expressed as ${\mathit{v}}_{{\varphi }_{\tau }}={[e^{j{\varphi }_{\tau ,1}},e^{j{\varphi }_{\tau 2}},\dots ,e^{j{\varphi }_{\tau ,N}}]}^T$.

In conclusion, considering the influence of the residual Doppler shift and the phase disturbance, the relationship between the real channel vector ${\mathit{h}}_{k,l}$ and the estimated channel vector ${\widehat{\mathit{h}}}_{k,l}$ can be expressed as:

$${\mathit{h}}_{k,l}={\widehat{\mathit{h}}}_{k,l}\odot {\mathit{v}}_{{\varphi }_{f_{d,k,l}^{rsd}}}\odot {\mathit{v}}_{{\varphi }_{\tau ,k,l}}=diag\left(diag\left({\widehat{\mathit{h}}}_{k,l}\right){\mathit{v}}_{{\varphi }_{f_{d,k,l}^{rsd}}}\right){\mathit{v}}_{{\varphi }_{{\tau }_{k,l}}},$$

(13)

where $\odot$ represents the Hadamard product. Let the channel phase of the $k\mathrm{th}$ user in the $l\mathrm{th}$ multicast group be ${\theta }_{k,l}={[{\theta }_{k,l,1},{\theta }_{k,l,2},\dots ,{\theta }_{k,l,N}]}^T$, which satisfies the uniform distribution between $0\sim 2\pi$. Then, the real channel phase with phase errors at instant $t_1$ is as follows:

$${\theta }_{k,l}(t_1)={\theta }_{k,l}(t_0)+{\varphi }_{f_{d,k,l}^{rsd}}+{\varphi }_{{\tau }_{k,l}}$$

(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 Figure 4 and Figure 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.

We adopt the similarity measurement function among multicast groups based on the $Ward$ connection method, i.e.,

$$d\left(i,j\right)=\sqrt{\frac{2n_i{n}_j}{n_i+n_j}}dist({\mathit{h}}_i^{eq},{\mathit{h}}_j^{eq}),$$

(15)

where $n_i$, $n_j$ represent the number of users of the group $i$ and the group $j$, respectively, ${\mathit{h}}_i^{eq}$, ${\mathit{h}}_j^{eq}$ represent the equivalent CSI of the group $i$ and the group $j$, respectively, and $dist({\mathit{h}}_i^{eq},{\mathit{h}}_j^{eq})$ represents the Euclidean distance between the vector ${\mathit{h}}_i^{eq}$ and the vector ${\mathit{h}}_i^{eq}$, i.e.,

$$dist({\mathit{h}}_i^{eq},{\mathit{h}}_j^{eq})=‖\frac{{\mathit{h}}_i^{eq}}{‖{\mathit{h}}_i^{eq}‖}-\frac{{\mathit{h}}_j^{eq}}{‖{\mathit{h}}_j^{eq}‖}‖$$

(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\mathrm{th}$ user in the $l\mathrm{th}$ multicast group can be expressed as:

$$\mathrm{{\rm E}}\left\{\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_{k,l}\right\}\approx \frac{\mathrm{{\rm E}}\left\{{\left|{\mathit{h}}_{k,l}^H{\mathit{f}}_l\right|}^2\right\}}{\mathrm{{\rm E}}\left\{{\displaystyle \sum _{j=1,j\ne l}^L{\left|{\mathit{h}}_{k,l}^H{\mathit{f}}_j\right|}^2}\right\}+{\sigma }^2},$$

(17)

$$R_{k,l}\approx B{\log }_2\left(1+\frac{\mathrm{{\rm E}}\left\{{\left|{\mathit{h}}_{k,l}^H{\mathit{f}}_l\right|}^2\right\}}{\mathrm{{\rm E}}\left\{{\displaystyle \sum _{j=1,j\ne l}^L{\left|{\mathit{h}}_{k,l}^H{\mathit{f}}_j\right|}^2}\right\}+{\sigma }^2}\right),$$

(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:

$$\mathrm{E}\mathrm{E}=\frac{{\displaystyle \sum _{l=1}^L\min \left({\left\{R_{k,l}\right\}}_{k=1}^{\left|U_l\right|}\right)}}{P_{total}}=\frac{B{\displaystyle \sum _{l=1}^L{\log }_2\left(1+\min {\left\{\frac{{\left|{\mathit{h}}_{k,l}^H{\mathit{f}}_l\right|}^2}{{\displaystyle \sum _{j=1,j\ne l}^L{\left|{\mathit{h}}_{k,l}^H{\mathit{f}}_j\right|}^2}+{\sigma }^2}\right\}}_{k=1}^{U_s}\right)}}{P_t+P_0},$$

(19)

where $P_{total}$ represents the total power consumption, $P_t={\displaystyle \sum _l^L\left|f_l{f}_l^H\right|}$ denotes the transmission power of the LEO satellite, and $P_0$ denotes the inherent power consumption of the communication system.

Let the Boolean variable ${\eta }_{k,l}\in \left\{0,1\right\}$ indicate whether the $k\mathrm{th}$ user in the $l\mathrm{th}$ multicast group is served, ${\eta }_{k,l}=1$ and ${\eta }_{k,l}=0$ indicate that the user can be served and not served, respectively, and $\mathit{\eta }=\left[{\mathit{\eta }}_1,{\mathit{\eta }}_2,\dots ,{\mathit{\eta }}_L\right]$, ${\mathit{\eta }}_l={\left[{\eta }_{1,l},{\eta }_{2,l},\dots ,{\eta }_{\left|U_l\right|,L}\right]}^T$. In conclusion, under the constraints of the transmission power and QoS, the problem of maximizing system EE can be modeled as:

$$Q_1: \underset{\mathit{\eta },{\mathit{f}}_l,{\mathrm{SINR}}_l^{\min }}{\max }\mathrm{E}\mathrm{E}=\frac{B{\displaystyle \sum _{l=1}^L{\log }_2\left(1+\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min }\right)}}{{\displaystyle \sum _l^L{\left|{\mathit{f}}_l{\mathit{f}}_l^H\right|}^2}+P_0},$$

(20)

$$s.t. C_1:{\eta }_{k,l}\in \left\{0,1\right\},\forall k,l,$$

(21)

$$C_2:\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_{k,l}\ge {\eta }_{k,l}\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min },\forall k,l,$$

(22)

$$C_3:\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min }\ge \mathrm{S}\mathrm{I}\mathrm{N}{R}_0,\forall l,$$

(23)

$$C_4:{\displaystyle \sum _{k=1}^{\left(U_l\right)}{\eta }_{k,l}}=U_s,\forall l,$$

(24)

$$C_5:{\displaystyle \sum _{l=1}^L{\left|{\mathit{f}}_l{\mathit{f}}_l^H\right|}^2}\le P_T,$$

(25)

where $\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min }$ represents the minimum SINR of the $l\mathrm{th}$ multicast group and constraint $C_3$ represents that $\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min }$ should be greater than the minimum SINR constraint $\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_0$. Constraint $C_4$ limits the number of scheduled users in each multicast group to $U_s$.

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_1$ 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_2$ and $C_5$ in the problem $Q_1$ involve the quadratic form of the variable ${\mathit{f}}_l$, therefore, $Q_1$ is the QCQP form problem. To handle this problem, we invoke the SDP algorithm, a new variable ${\mathit{W}}_l\triangleq {\mathit{f}}_l{\mathit{f}}_l^H$ is introduced, and the positive semidefinite matrix ${\mathit{W}}_l$ needs to meet the constraints of ${\mathit{W}}_l\underset{\_}{\succ }0$ and $rank\left({\mathit{W}}_l\right)=1$. Then, the problem $Q_1$ can be equivalent to:

$$Q_2: \underset{\mathit{\eta },{\mathit{W}}_l,\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min }}{\max }\mathrm{E}\mathrm{E}=\frac{B{\displaystyle \sum _{l=1}^L{\log }_2\left(1+\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min }\right)}}{{\displaystyle \sum _l^{L}Tr({\mathit{W}}_l)}+P_0},$$

(26)

$$s.t. C_1,C_2, C_3, C_4 in Q_1,$$

(27)

$$C_5:{\displaystyle \sum _l^{L}Tr({\mathit{W}}_l)}\le P_T,$$

(28)

$$C_6:{\mathit{W}}_l\underset{\_}{\succ }0, \forall l,$$

(29)

$$C_7:rank\left({\mathit{W}}_l\right)=1, \forall l,$$

(30)

Similarly, the SINR and the communication rate of the $k\mathrm{th}$ user in the $l\mathrm{th}$ multicast group can be equivalently converted to:

$$\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_{k,l}\approx \frac{\mathrm{{\rm E}}\left\{Tr\left({\mathit{H}}_{k,l}{\mathit{W}}_l\right)\right\}}{\mathrm{{\rm E}}\left\{{\displaystyle \sum _{j=1,j\ne l}^{L}Tr\left({\mathit{H}}_{k,l}{\mathit{W}}_j\right)}\right\}+{\sigma }^2},$$

(31)

$$\begin{array}{ll}{R}_{k,l}& =B{\log }_2\left(1+\frac{\mathrm{{\rm E}}\left\{Tr\left({\mathit{H}}_{k,l}{\mathit{W}}_l\right)\right\}}{\mathrm{{\rm E}}\left\{{\displaystyle \sum _{j=1,j\ne l}^{L}Tr\left({\mathit{H}}_{k,l}{\mathit{W}}_j\right)}\right\}+{\sigma }^2}\right)=B{\log }_2\left(1+\frac{Tr\left(\mathrm{{\rm E}}\left\{{\mathit{H}}_{k,l}{\mathit{W}}_l\right\}\right)}{{\displaystyle \sum _{j=1,j\ne l}^{L}Tr\left({\rm E}\left\{{\mathit{H}}_{k,l}{\mathit{W}}_j\right\}\right)}+{\sigma }^2}\right)\\ & =B{\log }_2\left(1+\frac{Tr\left({\overline{\mathit{H}}}_{k,l}{\mathit{W}}_l\right)}{{\displaystyle \sum _{j=1,j\ne l}^{L}Tr\left({\overline{\mathit{H}}}_{k,l}{\mathit{W}}_j\right)}+{\sigma }^2}\right)\end{array}$$

(32)

where ${\mathit{H}}_{k,l}\in C^{M\times M}$ is the instantaneous channel autocorrelation matrix of the $k\mathrm{th}$ user in the $l\mathrm{th}$ multicast group and ${\overline{\mathit{H}}}_{k,l}\in C^{M\times M}$ is the long-term channel autocorrelation matrix. The relationship between the two can be expressed as:

$${\overline{\mathit{H}}}_{k,l}=\mathrm{{\rm E}}\left\{{\mathit{H}}_{k,l}\right\}\triangleq \mathrm{{\rm E}}\left\{{\widehat{\mathit{h}}}_{k,l}{\widehat{\mathit{h}}}_{k,l}^H\right\}=diag\left({\widehat{\mathit{h}}}_{k,l}\right){\mathit{P}}_{f_{d,k,l}^{rsd}}{\mathit{Q}}_{{\tau }_{k,l}}diag\left({\widehat{\mathit{h}}}_{k,l}^H\right),$$

(33)

where ${\mathit{P}}_{f_{d,k,l}^{rsd}}$ and ${\mathit{Q}}_{{\tau }_{k,l}}$ can be expressed as follows:

$$\begin{array}{ll}{P}_{f_{d,k,l}^{rsd}}& ={\rm E}\left\{v_{{\varphi }_{f_{d,k,l}^{rsd}}}{v}_{{\varphi }_{f_{d,k,l}^{rsd}}}^H\right\}\\ & =\mathrm{{\rm E}}\left\{{\left[e^{j{\varphi }_{f_{d,k,l}^{rsd},1}},e^{j{\varphi }_{f_{d,k,l}^{rsd},2}},\dots ,e^{j{\varphi }_{f_{d,k.l}^{rsd},M}}\right]}^T\left[e^{-j{\varphi }_{f_{d,k,l}^{rsd},1}},e^{-j{\varphi }_{f_{d,k,l}^{rsd},2}},\dots ,e^{-j{\varphi }_{f_{d,k.l}^{rsd},M}}\right]\right\}\\ & =\mathrm{{\rm E}}\left\{\begin{array}{ccc}1& \cdots & e^{j{\varphi }_{f_{d,k,l}^{rsd},1}}{e}^{-j{\varphi }_{f_{d,k.l}^{rsd},M}}\\ ⋮& \ddots & ⋮\\ e^{j{\varphi }_{f_{d,k.l}^{rsd},M}}{e}^{-j{\varphi }_{f_{d,k,l}^{rsd},1}}& \cdots & 1\end{array}\right\}\\ & =\left\{\begin{array}{ccc}1& \cdots & {\rm E}\left\{e^{j{\varphi }_{f_{d,k,l}^{rsd},1}}{e}^{-j{\varphi }_{f_{d,k.l}^{rsd},M}}\right\}\\ ⋮& \ddots & ⋮\\ {\rm E}\left\{e^{j{\varphi }_{f_{d,k.l}^{rsd},M}}{e}^{-j{\varphi }_{f_{d,k,l}^{rsd},1}}\right\}& \cdots & 1\end{array}\right\}\end{array}$$

(34)

In (34), the diagonal elements of ${\mathit{P}}_{f_{d,k,l}^{rsd}}$ are all 1, and the elements in row $i$ and column $j$ on the non-diagonal are $\mathrm{{\rm E}}\left\{e^{j{\varphi }_{f_{d,k.l}^{rsd},i}}{e}^{-j{\varphi }_{f_{d,k,l}^{rsd},j}}\right\}=\mathrm{{\rm E}}\left\{e^{j{\varphi }_{f_{d,k.l}^{rsd},i}}\right\}\mathrm{{\rm E}}\left\{e^{-j{\varphi }_{f_{d,k,l}^{rsd},j}}\right\}$, according to ${\varphi }_{f_d^{rsd}}\sim N\left(0,{\sigma }_{{\varphi }_{f_d^{rsd}}}^2\right)$,

$$\begin{gathered} {\rm E}\left\{e^{j{\varphi }_{f_{d,k.l}^{rsd},i}}\right\}={{\displaystyle \int }}_{-\infty }^{\infty }{e}^{j{\varphi }_{f_{d,k.l}^{rsd},i}}\frac{1}{\sqrt{2\pi }{\sigma }_{{\varphi }_{f_{d,k,l}^{rsd}}}}{e}^{-\frac{{\varphi }_{f_{d,k,l}^{rsd}}^2}{2{\sigma }_{{\varphi }_{f_d^{rsd}}}^2}}d{\varphi }_{f_{d,k.l}^{rsd},i} \\ =e^{-\frac{{\sigma }_{f_{d,k,l}^{rsd}}^2}{2}}{{\displaystyle \int }}_{-\infty }^{\infty }\frac{1}{\sqrt{2\pi }{\sigma }_{{\varphi }_{f_{d,k,l}^{rsd}}}}{e}^{-\frac{{\left({\varphi }_{f_{d,k.l}^{rsd},i}-j{\sigma }_{{\varphi }_{f_{d,k,l}^{rsd}}}\right)}^2}{2{\sigma }_{{\varphi }_{f_d^{rsd}}}^2}}d{\varphi }_{f_{d,k.l}^{rsd},i} \end{gathered}$$

(35)

Similarly, $\mathrm{{\rm E}}\left\{e^{-j{\varphi }_{f_{d,k,l}^{rsd},j}}\right\}=e^{-\frac{{\sigma }_{{\varphi }_{f_d^{rsd}}}^2}{2}}$. Therefore, $\mathrm{{\rm E}}\left\{e^{j{\varphi }_{f_{d,k.l}^{rsd},i}}{e}^{-j{\varphi }_{f_{d,k,l}^{rsd},j}}\right\}=\mathrm{{\rm E}}\left\{e^{j{\varphi }_{f_{d,k.l}^{rsd},i}}\right\}{\rm E}\left\{e^{-j{\varphi }_{f_{d,k,l}^{rsd},j}}\right\}=e^{-{\sigma }_{{\varphi }_{f_d^{rsd}}}^2}$.

$$\begin{array}{ll}{\mathit{Q}}_{{\tau }_{k,l}}& =\mathit{{\rm E}}\left\{{\mathit{v}}_{{\varphi }_{{\tau }_{k,l}}}{\mathit{v}}_{{\varphi }_{t_{k,l}}}^H\right\}\\ & =\mathrm{{\rm E}}\left\{{\left[e^{j{\varphi }_{{\tau }_{k,l},1}},e^{j{\varphi }_{{\tau }_{k,l}2}},\dots ,e^{j{\varphi }_{{\tau }_{k,l},M}}\right]}^T\left[e^{-j{\varphi }_{{\tau }_{k,l},1}},e^{-j{\varphi }_{{\tau }_{k,l}2}},\dots ,e^{-j{\varphi }_{{\tau }_{k,l},M}}\right]\right\}\\ & =\mathrm{{\rm E}}\left\{\begin{array}{ccc}1& \cdots & e^{j{\varphi }_{{\tau }_{k,l},1}}{e}^{-j{\varphi }_{{\tau }_{k,l},M}}\\ ⋮& \ddots & ⋮\\ e^{j{\varphi }_{{\tau }_{k,l},M}}{e}^{-j{\varphi }_{{\tau }_{k,l},1}}& \cdots & 1\end{array}\right\}\\ & =\left\{\begin{array}{ccc}1& \cdots & {\rm E}\left\{e^{j{\varphi }_{{\tau }_{k,l},1}}{e}^{-j{\varphi }_{{\tau }_{k,l},M}}\right\}\\ ⋮& \ddots & ⋮\\ {\rm E}\left\{e^{j{\varphi }_{{\tau }_{k,l},M}}{e}^{-j{\varphi }_{{\tau }_{k,l},1}}\right\}& \cdots & 1\end{array}\right\}\end{array}$$

(36)

In (36), the diagonal elements of ${\mathit{Q}}_{{\tau }_{k,l}}$ are all 1, and the elements in row $i$ and column $j$ on the non-diagonal are $\mathrm{{\rm E}}\left\{e^{j{\varphi }_{{\tau }_{k,l},i}}{e}^{-j{\varphi }_{{\tau }_{k,l},j}}\right\}=\mathrm{{\rm E}}\left\{e^{j{\varphi }_{{\tau }_{k,l},i}}\right\}\mathrm{{\rm E}}\left\{e^{-j{\varphi }_{{\tau }_{k,l},j}}\right\}$, according to ${\varphi }_{\tau }\sim N\left(0,{\sigma }_{{\varphi }_{\tau }}^2\right)$,

$$\begin{array}{ll}\mathrm{{\rm E}}\left\{e^{j{\varphi }_{{\tau }_{k,l},i}}\right\}& ={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{e}^{j{\varphi }_{{\tau }_{k,l},i}}\frac{1}{\sqrt{2\pi }{\sigma }_{{\varphi }_{{\tau }_{k,l}}}^{}}}{e}^{-\frac{{\varphi }_{{}_{{\tau }_{k,l},i}}^2}{2{\sigma }_{{\varphi }_{{\tau }_{k,l}}}^2}}d{\varphi }_{{\tau }_{k,l},i}\\ & =e^{-\frac{{\sigma }_{{\varphi }_{{\tau }_{k,l}}}^2}{2}}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\frac{1}{\sqrt{2\pi }{\sigma }_{{\varphi }_{{\tau }_{k,l}}}}}{e}^{{}^{-\frac{{\left({\varphi }_{{\tau }_{k,l},i}-j{\sigma }_{{\varphi }_{{\tau }_{k,l}}}\right)}^2}{2{\sigma }_{{\varphi }_{{\tau }_{k,l}}}^2}}}d{\varphi }_{{\tau }_{k,l},i}\\ & =e^{-\frac{{\sigma }_{{\varphi }_{{\tau }_{k,l}}}^2}{2}}\end{array}$$

(37)

Similarly, $\mathrm{{\rm E}}\left\{e^{j{\varphi }_{{\tau }_{k,l},i}}\right\}=e^{-\frac{{\sigma }_{{\varphi }_{{\tau }_{k,l}}}^2}{2}}$. Therefore, $\mathrm{{\rm E}}\left\{e^{j{\varphi }_{{\tau }_{k,l},i}}{e}^{-j{\varphi }_{{\tau }_{k,l},i}}\right\}=\mathrm{{\rm E}}\left\{e^{j{\varphi }_{{\tau }_{k,l},i}}\right\}\mathrm{{\rm E}}\left\{e^{-j{\varphi }_{{\tau }_{k,l},i}}\right\}=e^{-{\sigma }_{{\varphi }_{{\tau }_{k,l}}}^2}$.

4.2. DC Programming

Since the constraint $C_1$ is a Boolean constraint and the constraint $C_2$ is a nonconvex constraint, the problem $Q_2$ is a nonconvex and nonsmooth combinatorial optimization problem. To handle this challenge, we can transform the problem $Q_2$ into a DC programming problem [33]. Therefore, the relaxation variable ${\zeta }_{k,l}$ is introduced as the lower bound of the SINR of the $k\mathrm{th}$ user in the $l\mathrm{th}$ multicast group, $\zeta =[{\zeta }_1,{\zeta }_2,\dots ,{\zeta }_L]$, ${\zeta }_l={[{\zeta }_{l,1},{\zeta }_{l,2},\dots ,{\zeta }_{l,\left|U_l\right|}]}^T$. The problem $Q_2$ can be equivalently converted to:

$$Q_3: \underset{\mathit{\eta },\mathit{W},\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min },\zeta }{\max }\mathrm{E}\mathrm{E}=\frac{B{\displaystyle \sum _{l=1}^L{\log }_2\left(1+\mathrm{S}\mathrm{I}\mathrm{N}{R}_l^{\min }\right)}}{{\displaystyle \sum _{l=1}^{L}Tr({\mathit{W}}_l)}+P_0},$$

(38)

$$s.t. C_1, C_3, C_4, C_5 ,C_6 ,C_7 in Q_1,$$

(39)

$$C_2:\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_{k,l}\ge {\zeta }_{k,l},\forall k,l,$$

(40)

$$C_8:{\zeta }_{k,l}\ge {\eta }_{k,l}\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min },\forall k,l,$$

(41)

where the constraint $C_2$ can be equivalently converted to:

$$C_2\Rightarrow 1+\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_{k,l}\ge 1+{\zeta }_{k,l},$$

(42)

To further express (42) in the form of DC programming, we introduce new function variables ${\mathrm{\Gamma }}_{k,l}\left(\mathit{W}\right)$ and ${\mathrm{{\rm I}}}_{k,l}\left(\mathit{W},{\zeta }_{k,l}\right)$:

$${\mathrm{\Gamma }}_{k,l}\left(\mathit{W}\right)={\sigma }^2+{\displaystyle \sum _{j=1.j\ne l}^{L}Tr\left({\mathit{H}}_{k,l}{\mathit{W}}_j\right)},$$

(43)

$${\mathrm{{\rm I}}}_{k,l}\left(\mathit{W},{\zeta }_{k,l}\right)=\frac{{\sigma }^2+{\displaystyle \sum _{j=1}^{L}Tr\left({\mathit{H}}_{k,l}{\mathit{W}}_j\right)}}{1+{\zeta }_{k,l}},$$

(44)

Therefore, the constraint $C_2$ can be rewritten as:

$$C_2\Rightarrow \mathrm{\Gamma }\left(\mathit{W}\right)-\mathrm{{\rm I}}\left(\mathit{W},{\zeta }_{k,l}\right)\le 0,$$

(45)

In (45), ${\mathrm{\Gamma }}_{k,l}\left({\mathit{W}}_l\right)$ is the affine function of $\mathit{W}$, ${{\rm I}}_{k,l}\left(\mathit{W},{\zeta }_{k,l}\right)$ is the concave function of $\mathit{W}$ and ${\zeta }_{k,l}$. The transformed constraint $C_2$ is a typical DC constraint.

Similarly, the constraint $C_8$ can be equivalently converted into the following DC form:

$$C_8\Rightarrow 4{\zeta }_{k,l}+{\left({\eta }_{k,l}-\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min }\right)}^2\ge {\left({\eta }_{k,l}+\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min }\right)}^2,$$

(46)

In (38), the objective function in the problem $Q_3$ 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_3$ into the following form:

$$Q_4:\underset{\mathit{\eta },\mathit{W},SIN{R}_l^{\min },\zeta }{\max }\mathrm{E}\mathrm{E}=2q{\left(B{\displaystyle \sum _{l=1}^L{\mathsf{\Upsilon }}_l(\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min })}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-q^2\left({\displaystyle \sum _{l=1}^{L}Tr({\mathit{W}}_l)}+P_0\right),$$

(47)

$$s.t. C_1, C_2,C_3, C_4, C_5 ,C_6 ,C_7,C_8 in Q_3,$$

(48)

where $q$ is the introduced auxiliary variable, and ${\mathsf{\Upsilon }}_l({\mathrm{SINR}}_l^{\min })$ is the introduced auxiliary function, which can be expressed as:

$${\mathsf{\Upsilon }}_l(\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min })={\log }_2\left(1+\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min }\right),$$

(49)

$$q=\frac{\sqrt{{\displaystyle \sum _{l=1}^L{\mathsf{\Upsilon }}_l({\mathrm{SINR}}_l^{\min })}}}{{\displaystyle \sum _{l=1}^{L}Tr({\mathit{W}}_l)}+P_0},$$

(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_1$ into $C_1\Rightarrow 0\le {\eta }_{k,l}\le 1, \forall k,l$. Meanwhile, to avoid the non-duality of the solution of ${\eta }_{k,l}$ caused by the relaxation, the penalty term $\mathrm{P}\left({\eta }_{k,l}\right)={\eta }_{k,l}\log {\eta }_{k,l}+\left(1-{\eta }_{k,l}\right)\log \left(1-{\eta }_{k,l}\right)$ is introduced into the objective function: let ${\lambda }_1>0$ be the penalty factor, and the problem $Q_4$ can be equivalently converted to:

$$Q_5: \underset{\mathit{\eta },\mathit{W},\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min },\zeta }{\max }\mathrm{E}\mathrm{E}=2q{\left(B{\displaystyle \sum _{l=1}^L{\mathsf{\Upsilon }}_l(\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min })}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-q^2\left({\displaystyle \sum _{l=1}^{L}Tr({\mathit{W}}_l)}+P_0\right)+{\lambda }_1{\displaystyle \sum _{l=1}^L\mathrm{P}\left({\eta }_{k,l}\right)},$$

(51)

$$s.t. C_2:{\mathrm{\Gamma }}_{k,l}\left(\mathit{W}\right)-{\mathrm{{\rm I}}}_{k,l}\left(\mathit{W},{\zeta }_{k,l}\right)\le 0,\forall k,l,$$

(52)

$$C_3, C_4, C_5 ,C_6 ,C_7 in Q_4,$$

(53)

$$C_8:4{\zeta }_{k,l}+{\left({\eta }_{k,l}-\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min }\right)}^2\ge {\left({\eta }_{k,l}+\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min }\right)}^2\forall k,l,$$

(54)

4.3. CCCP Algorithm

From (51), (52) and (54), it can be seen that the problem $Q_5$ 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_4$ of the problem $Q_5$ limits the selection of the initial point. To find a feasible initial point, we substitute the constraint $C_4$ into the objective function and set the penalty factor ${\lambda }_2>0$. Then, the problem $Q_5$ can be equivalently converted to:

$$\begin{array}{l}{Q}_6:\underset{\mathit{\eta },\mathit{W},\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min },\zeta }{\max }\mathrm{E}\mathrm{E}=2q{\left(B{\displaystyle \sum _{l=1}^L{\mathsf{\Upsilon }}_l(\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min })}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-q^2\left({\displaystyle \sum _{l=1}^{L}Tr({\mathit{W}}_l)}+P_0\right)\\ \hfill +{\lambda }_1{\displaystyle \sum _{l=1}^L\mathrm{P}\left({\eta }_{k,l}\right)}-{{\displaystyle \sum _l^L{\lambda }_2\left({\displaystyle \sum _k^{\left|U_l\right|}{\eta }_{k,l}-U_s}\right)}}^2\hfill \end{array}$$

(55)

$$s.t. C_2,C_3, C_5 ,C_6 ,C_7,C_8 in Q_5,$$

(56)

Convexification: Let ${\left({\eta }_{k,l},\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min },\mathit{W},\zeta \right)}^{\left(t-1\right)}$ be the estimated value of variables $\left({\eta }_{k,l},\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min },\mathit{W},\zeta \right)$ in iteration $t-1$ of the problem $Q_6$. In iteration $t$, for the convex part ${\lambda }_1{\displaystyle \sum _{l=1}^L\mathrm{P}\left({\eta }_{k,l}\right)}$, 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 $\mathrm{P}\left({\eta }_{k,l}\right)$ can be expressed as:

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

(57)

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

(58)

Similarly, in the constraint $C_2:{\mathrm{\Gamma }}_{k,l}\left(\mathit{W}\right)-{\mathrm{{\rm I}}}_{k,l}\left(\mathit{W},{\zeta }_{k,l}\right)\le 0$, we replace the concave function ${\mathrm{{\rm I}}}_{k,l}\left(\mathit{W},{\zeta }_{k,l}\right)$ with its first-order Taylor expansion, which is reflected in line three of Algorithm 1. The first-order Taylor expansion of ${\mathrm{{\rm I}}}_{k,l}\left(\mathit{W},{\zeta }_{k,l}\right)$ can be expressed as:

$${\mathrm{{\rm I}}}_{k,l}{\left(\mathit{W},{\zeta }_{k,l}\right)}^{te}={\mathrm{{\rm I}}}_{k,l}\left({\mathit{W}}^{\left(t-1\right)},{\zeta }_{k,l}^{\left(t-1\right)}\right)+{\nabla }^T{\mathrm{{\rm I}}}_{k,l}\left({\mathit{W}}^{\left(t-1\right)},{\zeta }_{k,l}^{\left(t-1\right)}\right)\left[\begin{array}{c}{\left\{{\mathit{W}}_l-{\mathit{W}}_l^{\left(t-1\right)}\right\}}_{l=1}^L\\ {\zeta }_{k,l}-{\zeta }_{k,l}^{\left(t-1\right)}\end{array}\right],$$

(59)

$$\nabla {\mathrm{{\rm I}}}_{k,l}\left({\mathit{W}}^{\left(t-1\right)},{\zeta }_{k,l}^{\left(t-1\right)}\right)={\left[{\left\{\frac{{\left({\overline{\mathit{H}}}_{k,l}\right)}^T}{1+{\zeta }_{k,l}^{\left(t-1\right)}}\right\}}_{l=1}^L,-\frac{{\sigma }^2+{\displaystyle \sum _{l=1}^{L}Tr\left({\overline{\mathit{H}}}_{k,l}{\mathit{W}}_l^{\left(t-1\right)}\right)}}{{\left(1+{\zeta }_{k,l}^{\left(t-1\right)}\right)}^2}\right]}^T,$$

(60)

In the constraint $C_8:4{\zeta }_{k,l}+{\left({\eta }_{k,l}-\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min }\right)}^2\ge {\left({\eta }_{k,l}+\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min }\right)}^2\forall k,l$, we replace the concave function ${\left({\eta }_{k,l}-\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min }\right)}^2$ with its first-order Taylor expansion, which is reflected in line 3 of Algorithm 1. The first-order Taylor expansion of ${\left({\eta }_{k,l}-\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min }\right)}^2$ can be expressed as:

$${\left({\eta }_{k,l}-\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min }\right)}^{2,te}={\left({\eta }_{k,l}^{\left(t-1\right)}-\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min ,\left(t-1\right)}\right)}^2+{\left[\begin{array}{c}2\left({\eta }_{k,l}^{\left(t-1\right)}-\mathrm{SIN}{\mathrm{R}}_l^{\min ,\left(t-1\right)}\right)\\ -2\left({\eta }_{k,l}^{\left(t-1\right)}-\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min ,\left(t-1\right)}\right)\end{array}\right]}^T\left[\begin{array}{c}{\eta }_{k,l}-{\eta }_{k,l}^{\left(t-1\right)}\\ \mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min }-\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min ,\left(t-1\right)}\end{array}\right],$$

(61)

Optimization: The optimization step is reflected in line nine of Algorithm 1. According to (57), (59) and (61), the problem $Q_6$ can be equivalently converted to

$$\begin{array}{ll}{Q}_7:\underset{\mathit{\eta },\mathit{W},\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min },\zeta }{\max }\mathrm{E}\mathrm{E}=& 2q{\left(B{\displaystyle \sum _{l=1}^L{\mathsf{\Upsilon }}_l(\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min })}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-q^2\left({\displaystyle \sum _{l=1}^{L}Tr({\mathit{W}}_l)}+P_0\right)\\ & +{\lambda }_1{\displaystyle \sum _{l=1}^L\mathrm{P}{\left({\eta }_{k,l}\right)}^{te}}-{{\displaystyle \sum _l^L{\lambda }_2\left({\displaystyle \sum _k^{\left|U_l\right|}{\eta }_{k,l}-U_s}\right)}}^2\hfill \end{array}$$

(62)

$$s.t. C_3, C_5 ,C_6 ,C_7 in Q_6,$$

(63)

$$C_2:{\mathrm{\Gamma }}_{k,l}\left(\mathit{W}\right)-{\mathrm{{\rm I}}}_{k,l}{\left(\mathit{W},{\zeta }_{k,l}\right)}^{te}\le 0,\forall k,l,$$

(64)

$$C_8:4{\zeta }_{k,l}+{\left({\eta }_{k,l}-\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min }\right)}^{2,te}\ge {\left({\eta }_{k,l}+\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min }\right)}^2\forall k,l,$$

(65)

For the problem $Q_7$, the variables ${\left({\eta }_{k,l},\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min },\mathit{W},\zeta \right)}^{\left(t+1\right)}$ 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.

• Initialize ${\eta }_{k,l}^{\left(0\right)}\approx 0$, $\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_0=1$;
• Find ${\mathit{W}}^{\left(0\right)}$, the following optimization problem are modeled:

  $$P_{FES}: \left\{{\mathit{W}}^{\left(0\right)}\right\}: \underset{\mathit{W}}{\min } {\displaystyle \sum _{l=1}^{L}Tr({\mathit{W}}_l)},$$

  (66)

  $$s.t. C_1:‖{\sigma }^2\dots Tr{({\overline{\mathit{H}}}_{k,l}{\mathit{W}}_j)}_{j\ne l}\dots ‖\le \frac{Tr({\overline{\mathit{H}}}_{k,l}{\mathit{W}}_l)}{{\eta }_{k,l}\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_0},\forall k,l,$$

  (67)

  $$C_2:{\displaystyle \sum _{l=1}^{L}Tr\left({\mathit{W}}_l\right)}\le P_T,$$

  (68)
• If $P_{FES}$ is feasible, proceed to the next step, otherwise, update ${\eta }_{k,l}^{\left(0\right)}=\delta {\eta }_{k,l}^{\left(0\right)},0<\delta <1$ and repeat step 2;
• Based on the ${\mathit{W}}^{\left(0\right)}$ obtained in step 2, calculate the $\mathrm{S}\mathrm{I}\mathrm{N}\mathrm{R}$ of each user, i.e., ${\mathrm{SINR}}_{k,l}^{\left(0\right)},\forall k,l$, and update ${\eta }_{k,l}^{\left(0\right)}$ according to ${\eta }_{k,l}^{\left(0\right)}=\min \left\{1,\frac{\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_{k,l}^{\left(0\right)}}{\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_0}\right\}$;
• Based on ${\eta }_{k,l}^{\left(0\right)}$ and ${\mathit{W}}^{\left(0\right)}$, calculate ${\zeta }^{\left(0\right)}$ and ${\left\{\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min ,}{}^{\left(0\right)}\right\}}_{l=1}^L$.

4.4. Penalty Iteration Algorithm

It should be noted that the SDP algorithm brings the nonconvex and nonsmooth constraint, i.e., $C_7:rank({\mathit{W}}_l)=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 ${\left\{{\mathit{W}}_l\right\}}_{l=1}^L$ meets the rank-one constraint. If so, the eigenvalue decomposition (EVD) algorithm is directly adopted to obtain the hybrid beamforming vectors ${\left\{{\mathit{f}}_l\right\}}_{l=1}^L$ according to ${\mathit{W}}_l={\mathit{f}}_l{f}_l^H$, and if the optimization solution ${\left\{{\mathit{W}}_l\right\}}_{l=1}^L$ 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 ${\left\{{\left\{{\mathit{w}}_{g,l}\right\}}_{g=1}^G\right\}}_{l=1}^L$ are generated based on the optimization solution ${\left\{{\mathit{W}}_l\right\}}_{l=1}^L$, 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({\mathit{W}}_l)=1$ is equivalent to $Tr({\mathit{W}}_l)-{\lambda }_{\max }\left({\mathit{W}}_l\right)=0$. Therefore, the nonsmooth method is adopted to transform the constraint $rank({\mathit{W}}_l)=1$ in the problem $Q_7$ into the following form:

$$C_7:Tr({\mathit{W}}_l)-{\lambda }_{\max }\left({\mathit{W}}_l\right)\le 0,$$

(69)

where ${\lambda }_{\max }\left({\mathit{W}}_l\right)$ is the function of solving the maximum eigenvalue. It should be noted that for any positive semidefinite matrix ${\mathit{W}}_l\underset{\_}{\succ }0$, the inequality $Tr({\mathit{W}}_l)-{\lambda }_{\max }\left({\mathit{W}}_l\right)\ge 0$ is always true, which means that the transformed constraint $C_7$ and $Tr({\mathit{W}}_l)-{\lambda }_{\max }\left({\mathit{W}}_l\right)=0$ are equivalent. Then, we can obtain that the matrix ${\mathit{W}}_l$ has only one non-zero eigenvalue and can be given by

$${\mathit{W}}_l={\lambda }_{\max }\left({\mathit{W}}_l\right){\mathit{w}}_{l,\max }{\mathit{w}}_{l,\max }^H,$$

(70)

where ${\mathit{w}}_{l,\max }$ is the corresponding unit eigenvector. Therefore, the problem $Q_7$ can be converted into

$$\begin{array}{ll}{Q}_8:\underset{\mathit{\eta },\mathit{W},\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min },\zeta }{\max }\mathrm{E}\mathrm{E}=& 2q{\left(B{\displaystyle \sum _{l=1}^L{\mathsf{\Upsilon }}_l(\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min })}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-q^2\left({\displaystyle \sum _{l=1}^{L}Tr({\mathit{W}}_l)}+P_0\right)\\ & +{\lambda }_1{\displaystyle \sum _{l=1}^L\mathrm{P}{\left({\eta }_{k,l}\right)}^{te}}-{{\displaystyle \sum _l^L{\lambda }_2\left({\displaystyle \sum _k^{\left|U_l\right|}{\eta }_{k,l}-U_s}\right)}}^2\end{array}$$

(71)

$$s.t. C_2, C_3, C_5 ,C_6 ,C_8 in Q_7,$$

(72)

$$C_7:Tr({\mathit{W}}_l)-{\lambda }_{\max }\left({\mathit{W}}_l\right)\le 0,$$

(73)

In the iterative calculation of the problem $Q_8$, based on the obtained ${\mathit{W}}_l$, if the value of $Tr({\mathit{W}}_l)-{\lambda }_{\max }\left({\mathit{W}}_l\right)$ is small enough, the matrix ${\mathit{W}}_l$ 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({\mathit{W}}_l)-{\lambda }_{\max }\left({\mathit{W}}_l\right)$ as small as possible, we adopt the penalty iteration algorithm and substitute the constraint $C_7$ into the objective function in the problem $Q_8$. Therefore, the problem $Q_8$ can be converted into

$$Q_9:\underset{\mathit{\eta },\mathit{W},\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min },\zeta }{\max }\mathrm{E}\mathrm{E}=F\left(\mathit{\eta },\mathit{W},\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min },\zeta \right)-{\lambda }_3{\displaystyle \sum _l^L\left(Tr\left({\mathit{W}}_l\right)-{\lambda }_{\max }\left({\mathit{W}}_l\right)\right)},$$

(74)

$$s.t. C_2, C_3, C_5 ,C_6 ,C_8 in Q_8$$

(75)

where $F\left(\mathit{\eta },\mathit{W},\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min },\zeta \right)=2q{\left(B{\displaystyle \sum _{l=1}^L{\mathsf{\Upsilon }}_l(\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min })}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-q^2\left({\displaystyle \sum _{l=1}^{L}Tr({\mathit{W}}_l)}+P_0\right)+{\lambda }_1{\displaystyle \sum _{l=1}^L\mathrm{P}{\left({\eta }_{k,l}\right)}^{te}}-{{\displaystyle \sum _l^L{\lambda }_2\left({\displaystyle \sum _k^{\left|U_l\right|}{\eta }_{k,l}-U_s}\right)}}^2$ and ${\lambda }_3>0$ is the penalty factor, which is generally larger enough to ensure that a smaller value of $Tr({\mathit{W}}_l)-{\lambda }_{\max }\left({\mathit{W}}_l\right)$ can be obtained.

According to (74), the iterative calculation of the problem $Q_9$ can maximize function $F\left(\mathit{\eta },\mathit{W},\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min },\zeta \right)$ and minimize function $Tr({\mathit{W}}_l)-{\lambda }_{\max }\left({\mathit{W}}_l\right)$. It should be noted that $Tr({\mathit{W}}_l)$ is an affine function and ${\lambda }_{\max }\left({\mathit{W}}_l\right)$ is nonsmooth, which can result in the nonsmoothness of the objective function in the problem $Q_9$. To handle the challenge, we replace ${\lambda }_{\max }\left({\mathit{W}}_l\right)$ with its first-order Taylor expansion. The subgradient of ${\lambda }_{\max }\left({\mathit{W}}_l\right)$ is $\frac{\partial {\lambda }_{\max }\left({\mathit{W}}_l\right)}{\partial {\mathit{W}}_l}={\mathit{w}}_{l,\max }{\mathit{w}}_{l,\max }^H$ and its first-order Taylor expansion can be expressed as follows, which is reflected in line eight of Algorithm 1.

$${\lambda }_{\max }\left({\mathit{W}}_l^{\left(t\right)}\right)\ge {\lambda }_{\max }\left({\mathit{W}}_l^{\left(t-1\right)}\right)+〈{\mathit{w}}_{l,\max }{\mathit{w}}_{l,\max }^H,{\mathit{W}}_l^{\left(t\right)}-{\mathit{W}}_l^{\left(t-1\right)}〉,$$

(76)

where $〈{\mathit{w}}_{l,\max }{\mathit{w}}_{l,\max }^H,{\mathit{W}}_l^{\left(t\right)}-{\mathit{W}}_l^{\left(t-1\right)}〉=Tr\left({\left({\mathit{w}}_{l,\max }{\mathit{w}}_{l,\max }^H\right)}^H\left({\mathit{W}}_l^{\left(t\right)}-{\mathit{W}}_l^{\left(t-1\right)}\right)\right)$.

We substitute (76) into the objective function in the problem $Q_9$ to replace ${\lambda }_{\max }\left({\mathit{W}}_l\right)$, and the problem $Q_9$ can be expressed as:

$$Q_{10}:\underset{\mathit{\eta },\mathit{W},\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min },\zeta }{\max }\mathrm{E}\mathrm{E}=F\left(\mathit{\eta },\mathit{W},\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min },\zeta \right)-{\lambda }_3{\displaystyle \sum _l^L\left(Tr\left({\mathit{W}}_l\right)-{\lambda }_{\max }\left({\mathit{W}}_l^{\left(t-1\right)}\right)+〈{\mathit{w}}_{l,\max }{\mathit{w}}_{l,\max }^H,{\mathit{W}}_l^{\left(t\right)}-{\mathit{W}}_l^{\left(t-1\right)}〉\right)},$$

(77)

$$s.t. C_2, C_3, C_5 ,C_6 ,C_8 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 ${\epsilon }_1$, penalty iteration algorithm iteration index $m$, thresholds ${\epsilon }_2$, penalty factor ${\lambda }_1$, ${\lambda }_2$, ${\lambda }_3$. 1. Initial: ${\left(\mathit{\eta },\mathit{W},\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min },\zeta \right)}^{\left(k=0\right)}$, $q^{\left(k=0\right)}$. 2. while $\left|{\mathrm{EE}}^{\left(k\right)}-{\mathrm{EE}}^{\left(k-1\right)}\right|\ge {\epsilon }_1$ 3.  Convexification step by (57), (59), (61). 4.  Calculation $q^{\left(k\right)}$, substitute $q^{\left(k\right)}$ into (77). 5.  Optimization step. 6.  Let ${\left(\mathit{\eta },\mathit{W},\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min },\zeta \right)}^{\left(m=0\right)}={\left(\mathit{\eta },\mathit{W},\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min },\zeta \right)}^{\left(k=0\right)}$. 7.  while $\left|{\left\{Tr\left({\mathit{W}}_l^{\left(m\right)}\right)-{\lambda }_{\max }\left({\mathit{W}}_l^{\left(m\right)}\right)\right\}}_{l=1}^L\right|\ge {\epsilon }_2$ 8.   Calculate the maximum eigenvalue ${\lambda }_{\max }\left({\mathit{W}}_l\right)$ of ${\mathit{W}}_l^{\left(m\right)}$ and the corresponding eigenvector ${\mathit{w}}_{l,\max }^{\left(m\right)}$. 9.   Using CVX toolbox, calculate the variables ${\left(\mathit{\eta },\mathit{W},\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min },\zeta \right)}_{opt}^{\left(m\right)}$ at the $m\mathrm{th}$ iteration according to (77). 10.    If ${\left\{{\mathit{W}}_l^{\left(m+1\right)}\right\}}_l^L\approx {\left\{{\mathit{W}}_l^{\left(m\right)}\right\}}_l^L$, then 11.     Update ${\lambda }_3=2{\lambda }_3$. 12.   else 13.     Update $m=m+1$. 14.   end 15.  end 16.  Update ${\left(\mathit{\eta },\mathit{W},\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min },\zeta \right)}^{\left(k+1\right)}={\left(\mathit{\eta },\mathit{W},\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min },\zeta \right)}^{\left(m\right)}$, $k=k+1$, ${\lambda }_1={\lambda }_1+1$, ${\lambda }_2={\lambda }_2+1$. 17. end Output:${\left(\mathit{\eta },\mathit{W},SIN{R}_l^{\min },\zeta \right)}_{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 $Q_{10}$ be ${\left(\mathit{\eta },\mathit{W},\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min },\zeta \right)}^{\left(k+1\right)}$ and $F_{EE}\left({\left(\mathit{\eta },\mathit{W},\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min },\zeta \right)}^{\left(k+1\right)}\right)$ at the $k\mathrm{th}$ iteration. Therefore, the convergence can be proved as follows:

$$\begin{array}{l}{F}_{EE}\left({\left(\mathit{\eta },\mathit{W},\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min },\zeta \right)}^{\left(k+1\right)}\right)=F\left({\left(\mathit{\eta },\mathit{W},\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min },\zeta \right)}^{\left(k+1\right)}\right)-{\lambda }_3{\displaystyle \sum _l^L\left(Tr\left({\mathit{W}}_l^{\left(k+1\right)}\right)-{\lambda }_{\max }\left({\mathit{W}}_l^{\left(k+1\right)}\right)\right)}\\ \ge F\left({\left(\mathit{\eta },\mathit{W},\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min },\zeta \right)}^{\left(k+1\right)}\right)-{\lambda }_3{\displaystyle \sum _l^L\left(Tr\left({\mathit{W}}_l^{\left(k+1\right)}\right)-{\lambda }_{\max }\left({\mathit{W}}_l^{\left(k\right)}\right)-〈{\mathit{w}}_{l,\max }{\mathit{w}}_{l,\max }^H,{\mathit{W}}_l^{\left(k+1\right)}-{\mathit{W}}_l^{\left(k\right)}〉\right)}\\ \stackrel{by\left(75\right)}{\ge }F\left({\left(\mathit{\eta },\mathit{W},\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min },\zeta \right)}^{\left(k\right)}\right)-{\lambda }_3{\displaystyle \sum _l^L\left(Tr\left({\mathit{W}}_l^{\left(k\right)}\right)-{\lambda }_{\max }\left({\mathit{W}}_l^{\left(k\right)}\right)\right)}\\ =F_{EE}\left({\left(\mathit{\eta },\mathit{W},\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min },\zeta \right)}^{\left(k\right)}\right)\end{array}$$

(79)

The convergence can be proved according to (79). Therefore, after initializing the values of ${\left(\mathit{\eta },\mathit{W},\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_l^{\min },\zeta \right)}^{\left(k=0\right)}$, ${\lambda }_1^{\left(k=0\right)}$, ${\lambda }_2^{\left(k=0\right)}$ and ${\lambda }_3^{\left(k=0\right)}$, 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\left(L{K}^{2}N\right)$. 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{\displaystyle \sum _{l=1}^L\left|U_l\right|}+2L$, the number of convex constraints is ${\displaystyle \sum _{l=1}^L\left|U_l\right|}$ and the number of linear constraints is ${\displaystyle \sum _{l=1}^L\left|U_l\right|}+3L$. Let the number of iterations in the penalty iteration algorithm and the CCCP algorithm be $I_p$ and $I_c$, respectively. In conclusion, the overall complexity of the proposed algorithm is $O\left(L{K}^{2}N+I_p{I}_c\left(2{\displaystyle \sum _{l=1}^L\left|U_l\right|}+2L\right)\left(2{\displaystyle \sum _{l=1}^L\left|U_l\right|}+3L\right)\right)$.

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, $U_{l,p}$ users in the $l\mathrm{th}$ multicast group are preselected as the user representatives, where $U_s<U_{l,p}\le \left|U_l\right|$. 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.

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., $\left|{\mathit{h}}_{k,l}^H{\mathit{f}}_l\right|$, the beamforming vector ${\mathit{f}}_l$ of the $l\mathrm{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 first step: according to the orthogonal criterion [11], a user is preselected for each multicast group in turn, which is reflected in line three and line four of Step 1 in Algorithm 2;
• The second step: based on the users of each multicast group selected in the first step, linearly correlated users are selected for each multicast group, which is reflected in line two and line three of Step 2 in Algorithm 2.

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 $\mathrm{I}{\mathrm{d}}_{\left(1\right)}=Index\left\{\max \left(‖{\mathit{h}}_{k,l}‖\right),\forall k,l\right\}$, select the user with the largest channel gain, $\mathrm{I}{\mathrm{d}}_{\left(1\right)}$ is the index of the user. 2. while $l\le L,l\ne \mathrm{I}{\mathrm{d}}_{\left(1\right)}$ 3.   For all users in the $l\mathrm{th}$ multicast group, calculate $Z_{k,l}={\mathit{h}}_{k,l}\left(I_N-{\displaystyle {\sum }_{j=\mathrm{I}{\mathrm{d}}_{\left(1\right)}}^{\mathrm{I}{\mathrm{d}}_{\left(l\right)}}\frac{{\mathit{h}}_{\left(j\right)}^H{\mathit{h}}_{\left(j\right)}}{{‖{\mathit{h}}_{\left(j\right)}‖}_2^2}}\right)$ in turn. 4.  $\mathrm{I}{\mathrm{d}}_{\left(l\right)}=Index\left\{\max \left(‖Z_{k,l}‖\right),\forall k\in l\right\}$, the user with index $\mathrm{I}{\mathrm{d}}_{\left(l\right)}$ is the preselected orthogonal user of the $l\mathrm{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\mathrm{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., $C_{k,l}={\mathit{h}}_{k,l}^H\frac{{\mathit{h}}_{\mathrm{I}{\mathrm{d}}_{\left(l\right)}}{\mathit{h}}_{\mathrm{I}{\mathrm{d}}_{\left(l\right)}}^H}{{‖{\mathit{h}}_{\mathrm{I}{\mathrm{d}}_{\left(l\right)}}^H‖}_2^2}$. 3.   Based on the $C_{k,l}$ of users in each multicast group, select top $\left(U_{l,p}-1\right)$ 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 ${\mathit{F}}_{BB}$ and analog beamforming matrix ${\mathit{F}}_{RF}$ in a hybrid beamformer. After obtaining $\mathit{W}$, we need to further solve ${\mathit{F}}_{BB}$ and ${\mathit{F}}_{RF}$. The solution method of ${\mathit{F}}_{BB}$ and ${\mathit{F}}_{RF}$ can be divided into two steps:

• The first step: we adopt the EVD algorithm to solve the hybrid beamforming matrix $\mathit{F}$ from $\mathit{W}$.
• The second step: we propose the MM-AltOpt algorithm to obtain ${\mathit{F}}_{BB}$ and ${\mathit{F}}_{RF}$.

7.1. Solution of The Hybrid Beamforming Matrix

Before calculating ${\mathit{F}}_{BB}$ and ${\mathit{F}}_{RF}$, it is necessary to obtain the hybrid beamforming matrix $\mathit{F}$. For the solution of $\mathit{F}$, we can adopt the EVD algorithm based on the previously obtained optimization variable $\mathit{W}$. According to the relationship between $\mathit{W}$ and $\mathit{F}$, i.e., ${\mathit{W}}_l\triangleq {\mathit{f}}_l{\mathit{f}}_l^H$, the solution of $\mathit{F}$ can be modeled as follows:

$$\underset{f_l}{\min }{‖{\mathit{W}}_l-{\mathit{f}}_l{\mathit{f}}_l^H‖}_F^2,$$

(80)

where the hybrid beamforming vector ${\mathit{f}}_l$ can be given by

$${\mathit{f}}_l=\sqrt{v_l}{\mathit{u}}_l,\forall l,$$

(81)

where $v_l$ is the maximum eigenvalue of the matrix ${\mathit{W}}_l$ and ${\mathit{u}}_l$ is the maximum eigenvector of the matrix ${\mathit{W}}_l$.

7.2. MM-AltOpt Algorithm: Solution of ${\mathit{F}}_{BB}$ and ${\mathit{F}}_{RF}$

According to $\mathit{F}={\mathit{F}}_{RF}{\mathit{F}}_{BB}$ and $\left|{\left({\mathit{F}}_{RF}\right)}_{i,j}\right|=1$, the solution of ${\mathit{F}}_{BB}$ and ${\mathit{F}}_{BB}$ can be modeled as a joint optimization problem with the power and constant modulus constraints, as follows:

$$P_1:\underset{{\mathit{F}}_{BB},{\mathit{F}}_{RF}}{\min } {‖{\mathit{F}}^{opt}-{\mathit{F}}_{RF}{\mathit{F}}_{BB}‖}_F^2,$$

(82)

$$s.t. C_1:{\mathit{F}}_{RF}\in \mathbb{F},$$

(83)

$$C_2:{‖{\mathit{F}}_{RF}{\mathit{F}}_{BB}‖}_F^2=P_T,$$

(84)

where $\mathbb{F}=\left\{{\mathit{F}}_{RF}\in C^{N\times L}|{\left|{\mathit{F}}_{RF}{}_{\left(i,j\right)}\right|}^2=1,1\le i\le N,1\le j\le L\right\}$ represents the unit modulus constraint, which is determined by the phase shifter in UPA, and the constraint $C_2$ represents the power constraint.

It is worth noting that the problem $P_1$ is a matrix decomposition problem with the constant modulus constraint and the equality constraint. The objective function is a nonconvex function of variables ${\mathit{F}}_{BB}$ and ${\mathit{F}}_{RF}$, and the constraints $C_1$ and $C_2$ 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 ${\mathit{F}}_{BB}$ to meet the power constraint. Then, for the solution of the analog beamforming matrix ${\mathit{F}}_{RF}$ with the unit modulus constraint, the MM algorithm is adopted [33].

7.2.1. 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 ${\mathit{F}}_{RF}$ based on the digital beamforming matrix $F_{BB}$. Thus, the problem $P_1$ can be expressed as:

$$P_2:\underset{{\mathit{F}}_{RF}}{\min } {‖\mathit{F}-{\mathit{F}}_{RF}{\mathit{F}}_{BB}^{\left(n\right)}‖}_F^2,$$

(85)

$$s.t. C_1:{\mathit{F}}_{RF}\in \mathbb{F},$$

(86)

where ${\mathit{F}}_{BB}^{\left(n\right)}$ represents the estimated value of the digital beamforming matrix ${\mathit{F}}_{BB}$ at the $n\mathrm{th}$ iteration. Due to the unit module constraint of elements in ${\mathit{F}}_{RF}$, the problem $P_2$ 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 $\mathit{F}$ by rows. According to the equivalence of the F-norm and $L_2$-norm of the vector, the problem $P_2$ can be rewritten as:

$$P_3: \underset{{\mathit{F}}_{RF}}{\min } {\displaystyle \sum _{i=1}^N{\mathit{F}}_i^H{\mathit{F}}_i}-2\Re \left({\mathit{F}}_i^H{\mathit{F}}_{BB}^{\left(n\right)}{\mathit{F}}_{RF,i}\right)+{\mathit{F}}_{RF,i}^H{\mathit{F}}_{BB}^{\left(n\right)}{\mathit{F}}_{BB}^{\left(n\right)}{}^H{\mathit{F}}_{RF,i},$$

(87)

$$s.t. C_1:{\mathit{F}}_{RF}\in \mathbb{F},$$

(88)

where ${\mathit{F}}_i^H$ represents the $i\mathrm{th}$ row vector of the matrix $\mathit{F}$, and ${\mathit{F}}_{RF,i}^H$ represents the $i\mathrm{th}$ row vector of the matrix ${\mathit{F}}_{RF}$.

It should be noted that the third term ${\mathit{F}}_{RF,i}^H{\mathit{F}}_{BB}^{\left(n\right)}{\mathit{F}}_{BB}^{\left(n\right)}{}^H{\mathit{F}}_{RF,i}$ in (86) is a convex function term, which needs further conversion. According to the first-order Taylor expansion, ${\mathit{F}}_{RF,i}^H{\mathit{F}}_{BB}^{\left(n\right)}{\mathit{F}}_{BB}^{\left(n\right)}{}^H{\mathit{F}}_{RF,i}$ can be converted into:

$$\begin{array}{ll}{\mathit{F}}_{RF,i}^H{\mathit{F}}_{BB}^{\left(n\right)}{\mathit{F}}_{BB}^{\left(n\right)}{}^H{\mathit{F}}_{RF,i}& ={\mathit{F}}_{RF,i}^{\left(q\right)H}{\mathit{F}}_{BB}^{\left(n\right)}{\mathit{F}}_{BB}^{\left(n\right)}{}^H{\mathit{F}}_{RF,i}^{\left(q\right)}+2\Re \left({\mathit{F}}_{RF,i}^{\left(q\right)H}{\mathit{F}}_{BB}^{\left(n\right)}{\mathit{F}}_{BB}^{\left(n\right)}{}^H\left({\mathit{F}}_{RF,i}-{\mathit{F}}_{RF,i}^{\left(q\right)}\right)\right)\\ & +{\left({\mathit{F}}_{RF,i}-{\mathit{F}}_{RF,i}^{\left(q\right)}\right)}^H{\mathit{F}}_{BB}^{\left(n\right)}{\mathit{F}}_{BB}^{\left(n\right)}{}^H\left({\mathit{F}}_{RF,i}-{\mathit{F}}_{RF,i}^{\left(q\right)}\right)\end{array}$$

(89)

where ${\mathit{F}}_{RF,i}^{\left(q\right)}$ represents the estimated value of ${\mathit{F}}_{RF,i}^{}$ at the $q\mathrm{th}$ iteration. According to the MM algorithm, the surrogate of ${\mathit{F}}_{RF,i}^H{\mathit{F}}_{BB}^{\left(n\right)}{\mathit{F}}_{BB}^{\left(n\right)}{}^H{\mathit{F}}_{RF,i}$ can be expressed as follows:

$$\begin{array}{ll}{\mathit{F}}_{RF,i}^H{\mathit{F}}_{BB}^{\left(n\right)}{\mathit{F}}_{BB}^{\left(n\right)}{}^H{\mathit{F}}_{RF,i}& \le {\mathit{F}}_{RF,i}^{\left(q\right)H}{\mathit{F}}_{BB}^{\left(n\right)}{\mathit{F}}_{BB}^{\left(n\right)}{}^H{\mathit{F}}_{RF,i}^{\left(q\right)}+2\Re \left({\mathit{F}}_{RF,i}^{\left(q\right)H}{\mathit{F}}_{BB}^{\left(n\right)}{\mathit{F}}_{BB}^{\left(n\right)}{}^H\left({\mathit{F}}_{RF,i}-{\mathit{F}}_{RF,i}^{\left(q\right)}\right)\right)\\ & +{\left({\mathit{F}}_{RF,i}-{\mathit{F}}_{RF,i}^{\left(q\right)}\right)}^H{\mathit{X}}^{\left(n\right)}\left({\mathit{F}}_{RF,i}-{\mathit{F}}_{RF,i}^{\left(q\right)}\right)\end{array}$$

(90)

where ${\mathit{X}}^{\left(n\right)}$ is a positive semidefinite matrix and satisfies the constraint ${\mathit{X}}^{\left(n\right)}\underset{\_}{\succ }{\mathit{F}}_{BB}^{\left(n\right)}{\mathit{F}}_{BB}^{\left(n\right)}{}^H$; here, we let ${\mathit{X}}^{\left(n\right)}={\lambda }_{\max }\left({\mathit{F}}_{BB}^{\left(n\right)}{\mathit{F}}_{BB}^{\left(n\right)}{}^H\right)\mathit{I}$ and ${\lambda }_{\max }\left({\mathit{F}}_{BB}^{\left(n\right)}{\mathit{F}}_{BB}^{\left(n\right)}{}^H\right)$ represents the maximum eigenvalue of the matrix ${\mathit{F}}_{BB}^{\left(n\right)}{\mathit{F}}_{BB}^{\left(n\right)}{}^H$. In conclusion, (89) can be further expressed as:

$$\begin{array}{ll}{\mathit{F}}_{RF,i}^H{\mathit{F}}_{BB}^{\left(n\right)}{\mathit{F}}_{BB}^{\left(n\right)}{}^H{\mathit{F}}_{RF,i}& \le {\lambda }_{\max }\left({\mathit{F}}_{BB}^{\left(n\right)}{\mathit{F}}_{BB}^{\left(n\right)}{}^H\right){\mathit{F}}_{RF,i}^H{\mathit{F}}_{RF,i}+2\Re \left({\mathit{F}}_{RF,i}^H\left({\mathit{F}}_{BB}^{\left(n\right)}{\mathit{F}}_{BB}^{\left(n\right)}{}^H-{\lambda }_{\max }\left({\mathit{F}}_{BB}^{\left(n\right)}{\mathit{F}}_{BB}^{\left(n\right)}{}^H\right)\mathit{I}\right){\mathit{F}}_{RF,i}\right)\\ & +{\mathit{F}}_{RF,i}^{\left(q\right)H}\left({\lambda }_{\max }\left({\mathit{F}}_{BB}^{\left(n\right)}{\mathit{F}}_{BB}^{\left(n\right)}{}^H\right)\mathit{I}-{\mathit{F}}_{BB}^{\left(n\right)}{\mathit{F}}_{BB}^{\left(n\right)}{}^H\right){\mathit{F}}_{RF,i}^{}\end{array}$$

(91)

According to (90), the surrogate function of the objective function in the problem $P_3$ can be expressed as:

$$\begin{gathered} P_3:\underset{F_{RF}}{min} {\displaystyle \sum }_{i=1}^N{F}_i^H{F}_i-2\Re \left(F_i^H{F}_{BB}^{\left(n\right)}{F}_{RF,i}\right)+F_{RF,i}^H{F}_{BB}^{\left(n\right)}{F}_{BB}^{\left(n\right)}{}^H{F}_{RF,i}\Rightarrow \\ {P}_4:\underset{F_{RF}}{min} {\displaystyle \sum }_{i=1}^N{F}_i^H{F}_i-2\Re \left(F_i^H{F}_{BB}^{\left(n\right)}{F}_{RF,i}\right)+\lambda {\left(F_{BB}^{\left(n\right)}{F}_{BB}^{\left(n\right)}{}^H\right)}_{RF,i}^H{{}_{RF,i}}_{max} \\ +2\Re \left(F_{RF,i}^H\left(F_{BB}^{\left(n\right)}{F}_{BB}^{\left(n\right)}{}^H-\lambda {\left(F_{BB}^{\left(n\right)}{F}_{BB}^{\left(n\right)}{}^H\right)}_{max}{()}_{RF,i}\right){()}_{RF,i}^{\left(q\right)H}\left(\lambda {\left(F_{BB}^{\left(n\right)}{F}_{BB}^{\left(n\right)}{}^H\right)}_{BB}^{\left(n\right)}{{}_{BB}^{\left(n\right)}}^H{}_{max}{()}_{RF,i}\right)\right), \end{gathered}$$

(92)

It is worth noting that the first and third terms of the objective function in the problem $P_4$ are constant terms, and the last term is independent of the variable ${\mathit{F}}_{RF,i}^H$. After ignoring the above three items, the problem $P_4$ can be converted to the following projection problem:

$$P_5: \underset{{\mathit{F}}_{RF}}{\min } {\displaystyle \sum _{i=1}^N{‖{\mathit{F}}_{RF,i}-{\mathit{c}}_i^{\left(q\right)}‖}_2^2},$$

(93)

$$s.t. C_1:{\mathit{F}}_{RF}\in \mathbb{F},$$

(94)

where ${\mathit{c}}_i^{\left(q\right)}={\mathit{F}}_{BB}^{\left(n\right)}{\mathit{F}}_i-\left({\mathit{F}}_{BB}^{\left(n\right)}{\mathit{F}}_{BB}^{\left(n\right)}{}^H-{\lambda }_{\max }\left({\mathit{F}}_{BB}^{\left(n\right)}{\mathit{F}}_{BB}^{\left(n\right)}{}^H\right)\mathit{I}\right){\mathit{F}}_{RF,i}^{\left(q\right)H}$.

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:

$${\mathit{F}}_{RF,i}^{}=e^{j\arg \left(c_i^{\left(q\right)}\right)},\forall i,$$

(95)

$${\mathit{F}}_{RF}^{}=e^{-j\arg \left(C^{\left(q\right)T}\right)},$$

(96)

where $C^{\left(q\right)}={\mathit{F}}_{BB}^{\left(n\right)}{\mathit{F}}^H-\left({\mathit{F}}_{BB}^{\left(n\right)}{\mathit{F}}_{BB}^{\left(n\right)}{}^H-{\lambda }_{\max }\left({\mathit{F}}_{BB}^{\left(n\right)}{\mathit{F}}_{BB}^{\left(n\right)}{}^H\right)\mathit{I}\right){\mathit{F}}_{RF}^{\left(q\right)H}$, which is reflected in line two of the inner algorithm in Algorithm 3.

7.2.2. Solution of The Digital Beamforming Matrix

Based on the analog beamforming matrix ${\mathit{F}}_{RF}$ obtained in the previous section, the solution problem of the digital beamforming matrix ${\mathit{F}}_{BB}$ can be modeled as follows:

$$P_6:\underset{{\mathit{F}}_{BB}}{\min } {‖\mathit{F}-{\mathit{F}}_{RF}^{\left(n\right)}{\mathit{F}}_{BB}‖}_F^2,$$

(97)

$$s.t. {‖{\mathit{F}}_{RF}^{\left(n\right)}{\mathit{F}}_{BB}‖}_F^2=P_T,$$

(98)

where ${\mathit{F}}_{RF}^{\left(n\right)}$ represents the estimated value of the analog beamforming matrix ${\mathit{F}}_{RF}$ at the $n\mathrm{th}$ iteration. Due to the quadratic form and the convex equality constraint in the problem $P_6$, the problem $P_6$ is a nonconvex QCQP form problem.

One of the ways to solve the problem $P_6$ is to relax the equality constraint into the inequality constraint, and then convert the problem $P_6$ 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 $\mathit{F}$ satisfies the power constraint, i.e., ${‖\mathit{F}‖}_F^2=P_T$, 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:

$${\mathit{F}}_{BB}={\left({\mathit{F}}_{RF}^{\left(n\right)H}{\mathit{F}}_{RF}^{\left(n\right)}\right)}^{-1}{\mathit{F}}_{RF}^{\left(n\right)H}\mathit{F},$$

(99)

In conclusion, the MM-AltOpt algorithm for solving ${\mathit{F}}_{BB}$ and ${\mathit{F}}_{RF}$ 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 $\mathit{F}$, initial: ${\mathit{F}}_{RF}^{\left(n=0\right)}\in \mathbb{F}$, iteration index $n=0$, threshold ${\epsilon }_3={10}^{-3}$, the solution of the objective function of the problem $P_1$ in the $n\mathrm{th}$ iteration is ${\delta }^{\left(n\right)}$. 1. while $\left|{\delta }^{\left(n\right)}-{\delta }^{\left(n-1\right)}\right|\ge {\epsilon }_3$ 2.   Based on ${\mathit{F}}_{RF}^{\left(n\right)}$, calculate ${\mathit{F}}_{BB}^{\left(n+1\right)}$ according to ${\mathit{F}}_{BB}={\left({\mathit{F}}_{RF}^{\left(n\right)H}{\mathit{F}}_{RF}^{\left(\left(n\right)\right)}\right)}^{-1}{\mathit{F}}_{RF}^{\left(n\right)H}\mathit{F}$. 3.   Based on ${\mathit{F}}_{BB}^{\left(n+1\right)}$, calculate ${\mathit{F}}_{RF}^{\left(n+1\right)}$ according to the inner algorithm. 4.   Set $n=n+1$. 5. end Output: ${\mathit{F}}_{RF}$,${\mathit{F}}_{BB}$, normalize ${\mathit{F}}_{BB}=\frac{\sqrt{P_T}}{{‖{\mathit{F}}_{{}_{RF}}{\mathit{F}}_{BB}‖}_{\mathit{F}}}{\mathit{F}}_{BB}$. Inner algorithm: Algorithm for solving the analog beamforming matrix. Input: Hybrid beamforming matrix $\mathit{F}$, ${\mathit{F}}_{BB}^{\left(n\right)}$, ${\mathit{F}}_{RF}^{\left(q=0\right)}\in \mathbb{F}$, iteration index $q=0$, threshold ${\epsilon }_4={10}^{-3}$, the solution of the objective function of the problem $P_5$ in the $q\mathrm{th}$ iteration is ${\delta }_1^{\left(q\right)}$. 1. while $\left|{\delta }^{\left(n\right)}-{\delta }^{\left(n-1\right)}\right|\ge {\epsilon }_4$ 2.     Calculate $C^{\left(q\right)}={\mathit{F}}_{BB}^{\left(n\right)}{\mathit{F}}^H-\left({\mathit{F}}_{BB}^{\left(n\right)}{\mathit{F}}_{BB}^{\left(n\right)}{}^H-{\lambda }_{\max }\left({\mathit{F}}_{BB}^{\left(n\right)}{\mathit{F}}_{BB}^{\left(n\right)}{}^H\right)\mathit{I}\right){\mathit{F}}_{RF}^{\left(q\right)H}$. 3.     Calculate ${\mathit{F}}_{RF}^{}=e^{-j\arg \left(C^{\left(q\right)T}\right)}$. 4.     Set $q=q+1$. 5. end Output: ${\mathit{F}}_{RF}^{\left(n\right)}$.

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 $\mathrm{S}\mathrm{I}\mathrm{N}{\mathrm{R}}_0=1$. To facilitate analysis, we assume that the CSI errors of different multicast groups are the same, which are expressed as ${\sigma }_{f_{d,k,l}^{rsd}}^2={\sigma }_{f_d^{rsd}}^2$ and ${\sigma }_{{\varphi }_{{\tau }_{k,l}}}^2={\sigma }_{{\varphi }_{\tau }}^2$. The value of $P_0$ can be calculated by [31]. In addition, the system parameters used in the numerical simulations are shown in

Table 1. Simulation parameters.

Parameters	Values	Parameters	Values
$N$	$8\times 8=64$	$\kappa$	$1.38\times {10}^{-23}\mathrm{J}·{\mathrm{K}}^{-1}$
$L$	7	$P_0$	$21.5 \mathrm{W}$
${\kappa }_{k,l}$	10	$T$	300 K
Bandwidth	50 MHz	$G_{leo}$	3 dB
Orbit altitude	1000 km	$G_{ut}$	3 dB
Beam radius	250 km	$f$	20 GHz
$L{P}_{at}$	0.017 dB	$K$	150

.

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., ${\sigma }_{f_d^{rsd}}^2=25$, ${\sigma }_{{\varphi }_{\tau }}^2=10$ and ${\sigma }_{f_d^{rsd}}^2=20$, ${\sigma }_{{\varphi }_{\tau }}^2=5$. In addition, we set $U_s=2$ and $P_T=50 \mathrm{W}$. Meanwhile, we set two different numbers of the preselected users, i.e., $U_{l,p}/U_s=2,U_{l,p}=4$ and $U_{l,p}/U_s=3,U_{l,p}=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 ${\sigma }_{f_d^{rsd}}^2=25$, ${\sigma }_{{\varphi }_{\tau }}^2=10$ and ${\sigma }_{f_d^{rsd}}^2=20$, and ${\sigma }_{{\varphi }_{\tau }}^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 8 compares the EE performance of the proposed algorithm and the traditional algorithm under different system parameters, versus different transmission power thresholds $P_T$. 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 $P_T$. Under the conditions of ${\sigma }_{f_d^{rsd}}^2=25$, ${\sigma }_{{\varphi }_{\tau }}^2=10$ and $P_T=15 \mathrm{W}$, when $U_{l,p}/U_s=2,U_{l,p}=4$ and $U_{l,p}/U_s=3,U_{l,p}=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.

Figure 9 indicates the change trend of the SE of the massive MIMO LEO satellite multigroup multicast communication system versus different transmission power thresholds $P_T$. 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 ${\sigma }_{f_d^{rsd}}^2=25$, ${\sigma }_{{\varphi }_{\tau }}^2=10$ and $P_T=30 \mathrm{W}$, when $U_{l,p}/U_s=2,U_{l,p}=4$ and $U_{l,p}/U_s=3,U_{l,p}=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.

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 11 shows the change trajectory of the system EE and SE versus the different $U_s$. It can be seen that with the increase in $U_s$, 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_s$, 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_s$, 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_s$, 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 $U_{l,p}/U_s=2,U_{l,p}=4$, when $U_s=2$, the EE performance gain is 21.05% higher than when $U_s=4$.

Figure 12 shows the performance comparison of different algorithms for solving ${\mathit{F}}_{BB}$ and ${\mathit{F}}_{RF}$. 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 $\mathit{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_T=10 \mathrm{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\left(I_{MM}\left(N^3+I_{Inner}2N{L}^2+NL\right)\right)$, $O\left(I_{PE}\left(N^3+L^3+NL\right)\right)$ and $O\left(I_{OMP}\left(L^{4}N+L^2+N^2{L}^2+2L^3\right)\right)$, respectively, where $I_{MM}$, $I_{Inner}$, $I_{PE}$ and $I_{OMP}$ 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.

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.