## 4 Joint User Scheduling and Hybrid Beamforming Design for Maximizing EE
### Penalty Iteration Algorithm  
It should be noted that the SDP algorithm brings the nonconvex and nonsmooth constraint, i.e., \(C_{7}:rank(\mathbf{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 \(\{\mathbf{W}_{l}\}_{l=1}^{L}\) meets the rank-one constraint. If so, the eigenvalue decomposition (EVD) algorithm is directly adopted to obtain the hybrid beamforming vectors \(\{f_{l}\}_{l=1}^{L}\) according to \(\mathbf{W}_{l}=f_{l}f_{l}^{H}\), and if the optimization solution \(\{\mathbf{W}_{l}\}_{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 \(\bigg{\{}\Big{\{}\boldsymbol{w}_{gJ}\Big{\}}_{g=1}^{G}\bigg{\}}_{l=1}^{L}\) are generated based on the optimization solution \(\{\mathbf{W}_{l}\}_{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(\mathbf{W}_{l})=1\) is equivalent to \(Tr(\mathbf{W}_{l})-\lambda_{\text{max}}(\mathbf{W}_{l})=0\). Therefore, the nonsmooth method is adopted to transform the constraint \(rank(\mathbf{W}_{l})=1\) in the problem \(Q_{7}\) into the following form:  
\[C_{7}:Tr(\mathbf{W}_{l})-\lambda_{\text{max}}(\mathbf{W}_{l})\leq 0, \tag{69}\]  
where \(\lambda_{\text{max}}(\mathbf{W}_{l})\) is the function of solving the maximum eigenvalue. It should be noted that for any positive semidefinite matrix \(\mathbf{W}_{l}\succeq 0\), the inequality \(Tr(\mathbf{W}_{l})-\lambda_{\text{max}}(\mathbf{W}_{l})\geq 0\) is always true, which means that the transformed constraint \(C_{7}\) and \(Tr(\mathbf{W}_{l})-\lambda_{\text{max}}(\mathbf{W}_{l})=0\) are equivalent. Then, we can obtain that the matrix \(\mathbf{W}_{l}\) has only one non-zero eigenvalue and can be given by  
\[\mathbf{W}_{l}=\lambda_{\text{max}}(\mathbf{W}_{l})\mathbf{w}_{l,\text{max}}\mathbf{w}_{l, \text{max}}^{H}, \tag{70}\]  
where \(\mathbf{w}_{l,\text{max}}\) is the corresponding unit eigenvector. Therefore, the problem \(Q_{7}\) can be converted into  
\[Q_{8}:\max_{\mathbf{\eta},\text{W},\text{SINR}_{l}^{\text{min}}, \zeta}\text{EE} = 2q\bigg{(}B\sum\limits_{l=1}^{L}Y_{l}(\text{SINR}_{l}^{\text{ min}})\bigg{)}^{\frac{1}{2}}-q^{2}\bigg{(}\sum\limits_{l=1}^{L}Tr(\mathbf{W}_{l})+P_{0} \bigg{)} \tag{71}\] \[+\lambda_{1}\sum\limits_{l=1}^{L}\text{P}(\eta_{k,l})^{te}-\sum \limits_{l}^{L}\lambda_{2}\bigg{(}\sum\limits_{k}^{|U_{l}|}\eta_{k,l}-U_{s} \bigg{)}^{2}\] \[s.t.\ C_{2},\ C_{3},\ C_{5},C_{6},C_{8}\ in\ Q_{7},\] (72) \[C_{7}:Tr(\mathbf{W}_{l})-\lambda_{\text{max}}(\mathbf{W}_{l})\leq 0, \tag{73}\]  
In the iterative calculation of the problem \(Q_{8}\), based on the obtained \(\mathbf{W}_{l}\), if the value of \(Tr(\mathbf{W}_{l})-\lambda_{\text{max}}(\mathbf{W}_{l})\) is small enough, the matrix \(\mathbf{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(\mathbf{W}_{l})-\lambda_{\text{max}}(\mathbf{W}_{l})\) 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}\).