Robust Energy-Efficient Hybrid Beamforming Design for Massive MIMO LEO Satellite Communication Systems

In low earth orbit (LEO) satellite communication systems, the limited energy supply capacity and the difficulty in obtaining channel state information (CSI) are the practical challenges. Motivated by this, we focus on the design of robust hybrid beamforming to maximize the energy efficiency of LEO satellite communication systems. Assuming that the LEO satellite transmitter adopts the massive multi-input multi-output (MIMO) technology, considering the CSI errors caused by propagation delay and Doppler shift, under the constraints of transmit power and quality of service (QoS), a robust energy-efficient hybrid beamforming scheme is proposed. Since that there are no explicit expressions for the ergodic user rate and the ergodic signal-to-interference-plus-noise ratio, the approximate values with closed-form expressions are adopted. Then, we invoke the semidefinite programming (SDP) algorithm to transform the nonconvex quadratic constrained quadratic programming (QCQP) problem equivalently, and an inner and outer nested iterative algorithm combining quadratic transformation fractional programming (QTFP) and concave convex process (CCCP) is utilized to transfer a nonconvex problem into a convex problem. Meanwhile, we adopt a penalty function algorithm to solve the rank-one constraint in semidefinite programming algorithm. Finally, we invoke the normal form distance minimization (NFDM) algorithm and the alternating optimization (AltOpt) algorithm to jointly solve the digital beamforming matrix and analog beamforming matrix in hybrid beamformer. Numerical results validate that our proposed robust approach significantly outperforms the conventional one.

can not meet the future high-speed rate communication requirement. Therefore, we comprehensively consider the factors of energy efficiency and system capacity, the hybrid beamforming architecture based on the full connection structure is a cost-effective choice [14]. ·Due to the limited energy supply capacity on the LEO satellites, the transmit power constraint of the transmitter should be considered. And the QoS of various terminals are different, the signal-tointerference-plus-noise ratio (SINR) constraint of each user should be taken into account. In conclusion, aiming at maximizing the EE, we focus on the robust energy-efficient hybrid beamforming design for LEO satellite communication systems. The major contributions of the current work are summarized as follows: ·We introduce Massive MIMO transmission technology into LEO satellite communication systems. And assuming that the satellite transmitter is equipped with a uniform plane array (UPA), we establish the downlink model of LEO satellite communication systems. · We establish the Massive MIMO channel model for LEO satellite communication systems by incorporating the LEO satellite signal propagation properties. Meanwhile, we analyze the effects of propagation delay and Doppler shift on the CSI. ·Due to the CSI errors, the ergodic user rate and ergodic SINR can not be expressed explicitly. To this end, we adopt the explicit tight approximations of the two ergodic expressions. · For the quadratic constrained quadratic programming (QCQP) optimization problem, to simplify the problem modeling, we adopt the semidefinite programming (SDP) algorithm to transform the optimization problem equivalently. · To solve the problem of sum-of-ratios fractional programming form and the nonconvex characteristic in the objective function, we propose an inner and outer nested iteration algorithm combining quadratic transformation fractional programming (QTFP) and concave convex process (CCCP) to transfer the objective function into a concave function. Then, we adopt the convex optimization algorithm to solve the problem. · To handle the rank-one constraint problem existing in the SDP algorithm, we adopt a penalty function algorithm. · To obtain the digital beamforming matrix and the analog beamforming matrix in the hybrid beamformer, we adopt the low-complexity alternating optimization (AltOpt) algorithm on the basis of the normal form distance minimization (NFDM) algorithm. · Simulation results show that the robust algorithm proposed in this paper has a fast convergence speed and a significant performance gain compared with the conventional non-robust algorithm.
·For LEO satellite communication systems, we introduce the Massive MIMO transmission technology, which can improve the spatial resolution, spectral efficiency and energy efficiency. Meanwhile, aiming at the influence of CSI errors, we propose a robust method, which can improve the reliability of the communication system. Therefor, the proposed method can be suitable for the high-frequency band, high bandwidth, high-speed rate and high mobility communication scenario. Note that The major notations adopted in the paper is listed in Table I for ease of reference.  TABLE I

II. SYSTEM MODEL
As depicted in Fig.1 , the spectral efficiency of the hybrid beamforming system is comparable to that of the digital beamforming system [15]. In this paper, we mainly investigate the hybrid beamforming design under the The data streams at the LEO satellite transmitter are first processed by the digital beamformer for baseband processing, and then mapped to the UPA through the RF links for analog beamforming. Let   where f represents the carrier frequency, k L represents the number of channel propagation paths of the th k user, , , represent the complex gain, Doppler shift, propagation delay, and UPA array response vector of path l , respectively. In the following, we mainly focus on three factors: Doppler shift, propagation delay, and array response.

A. DOPPLER SHIFT
The Doppler shift ,  [19]. In practical communication systems, the Doppler shift is usually estimated by the frequency offset estimation method, and then Doppler shift precompensation is implemented [20]. The frequency offset estimation is generally divided into two stages: coarse estimation and fine estimation. During the coarse estimation stage, the user terminals calculate the Doppler shift at each moment based on the LEO satellite ephemeris and their own position information. As depicted in Fig. 2, the calculation formula is defined as: where , c f represent the speed of light and carrier frequency, respectively, v represents the vector of satellite motion velocity, l  represents the angle between the satellite motion direction and the signal propagation path l , , R H represent the earth's radius and orbital altitude, respectively,  represents the elevation angle of the user associated with path l . The coarse estimation can address the problem of large Doppler shift caused by the relative motion between the LEO satellite and user terminals. During the fine estimation stage, the users further calculate the Doppler shift ,

C. ARRAY RESPONSE
The array response vector of the UPA in (2) can be expressed as [25]: , , where  represents Kronecker product, In the high mobility communication scene, according to the LEO satellite trajectory and the general law of terminal motion, we can calculate the spatial angle information of the next moment at insatnt t to avoid the expiration of the information.
In conclusion, the channel vector ( , ) k t f h can be further expressed as: where ( , ) k g t f represents the downlink channel gain of the th k user, which can be expressed as The statistical characteristic of ( , ) k g t f depends on the propagation environment, for the LEO satellite communication scenario, the link power loss mainly includes the free space path loss and the atmospheric absorption loss. The free space path loss f s LP can be given by 10 10 10 where d represents the the traveled distance, c represents the speed of light. The atmospheric absorption loss at LP can refer to the reference [44].
Because of line of sight (LoS) transmission in LEO satellite communication systems, ( , ) k g t f can be modeled by Rician fading [26]. Let the Rician factor be k k and the power Considering the influence of residual Doppler shift and residual propagation delay in the channel estimation, let the estimated channel vector beˆk h and the actual channel vector be k h . Then, we model the actual channel vector k h as:( where e represents Hadamard product, represents the channel error vector caused by residual propagation delay, which can be expressed as: Then the actual channel phase of the th k user at 1 t can be expressed as: represent the channel phase components with the elements independently and uniformly distributed between 0 and 2 .

IV. PROBLEM FORMULATION
Aiming at maximizing the EE, we investigate the robust energy-efficient hybrid beamforming design for LEO satellite communication systems. The EE can be denoted as the ratio of system rate to total power consumption.

A. SYSTEM RATE
The SINR of the th k user in the system can be calculated according to (1), i.e., According to the analysis of the channel model in the previous section, the actual CSI of the LEO satellite communication system is a constantly changing random process, and it is not feasible to obtain the instantaneous and accurate CSI. To this end, we adopt the statistics averaging method for the system rate modelling [27], and the ergodic capacity of the LEO satellite communication system can be expressed as: where B represents the channel bandwidth, and (

B. SYSTEM POWER CONSUMPTION
The power consumption of the LEO satellite communication system includes transmit power consumption and circuit hardware power consumption. The total power consumption can be expressed as: This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication.
where the first term is the power of the th k user radiated by the transmitting antennas, which should meet the transmit power constraint of the LEO satellite communication system, i.e., The second term 0 P represents the power consumption of circuit hardware, such as the transmitting antennas in the LEO satellite communication system. In this paper, the LEO satellite transmitter adopts a hybrid beamforming architecture based on full connection, and each user terminal is configured with a single antenna. Thus, 0 P can be expressed as [28]: represent the power consumption of phase shifter, local oscillator, precoder, digital to analog converter, mixer, low-pass filter, and baseband amplifier, respectively.

C. PROBLEM MODELING
In conclusion, the system EE can be defined as: Considering the constraints of transmit power and QoS, and taking EE maximization as the optimization goal, the optimization problem can be modeled as: where k r represents the SINR constraint threshold of the th k user.

V. ROBUST HYBRID BEAMFORMING DESIGN
It is worth noting that the ergodic user rate and ergodic SINR would not allow explicit expressions, to handle this challenge, we invoke the explicit tight approximations of them. To simplify problem 1 Q , we adopt SDP algorithm to transform the optimization problem equivalently. We can observe that optimization problem is a general nonconvex sum-of-ratios fractional programming problem, to this end, we propose an inner and outer nested iterative algorithm combining QTFP and CCCP. In addition, it is worth noting that there is a rankone constraint in SDP algorithm, which can be handled by adopting the penalty function algorithm. Finally, we adopt low complexity NFDM algorithm and AltOpt algorithm to jointly solve digital beamforming matrix and analog beamforming matrix.

A. APPROXIMATED ERGODIC USER RATE AND ERGODIC SINR
We can observe that both the ergodic user rate and the ergodic SINR in (15) and (14) do not admit explicit expressions, to handle this problem, some researches use Monte Carlo method [29] for statistical simulation, which might be not practical due to the high computational complexity and huge demand for storage memory. To handle this challenge, we invoke the approximation k R of ergodic user rate and the approximation k SINR of ergodic SINR as follows [30]: It should be noted that the approximations in (26) and (27) are both very tight and have been theoretically analyzed and numerically verified in previous works [30].

B. SDP ALGORITHM
We can observe that the objective function and constraints (24) , (25) in problem 1 Q involve the second power of k b , thus, problem 1 Q is a typical nonconvex quadratic constrained quadratic programming (QCQP) problem [31]. To handle this challenge, we adopt the SDP algorithm [32], and then transform the optimization variables   1 This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication.
where k  Q represents the expectation of the autocorrelation matrix of the channel error vector It can be observed that the diagonal elements in Besides, k f p follows the distribution of Jacks , and its autocorrelation matrix k f P can be expressed as: . . (35), (36), (37), (38) s t

C. QTFP ALGORITHM
We can observe that 3 Q is a sum-of-ratios fractional programming problem, to this end, we adopt the quadratic transformation method. According to the quadratic transformation theory [33], problem 3 Q is equivalent to . . (35), (36), (37)

D. CCCP ALGORITHM
We can observe from (41) and (42)  are both concave functions of W , so 3 Q is a difference of convex (DC) programming [34]. To handle this challenge, we invoke the CCCP algorithm to address this DC programming. The CCCP algorithm is a monotonically decreasing global optimization method, which can be used to solve nonconvex optimization problems [35]. The basic idea is to replace (42) After transformation, the objective function in problem 5 Q is transformed into a concave function, which can be solved by convex optimization method. The basic idea is to initialize the feasible solution   0   W , use the auxiliary variable y in the iteration, and bring it into 5 Q for iterative solution, until the convergence condition or the maximum number of iterations is reached.

1) ALGORITHM APPLICATION ANALYSIS
It is worth noting that there is a rank-one constraint in problem 5 Q , to solve this problem, some references adopt semidefinite relaxation (SDR) algorithm [36], which directly removes the rank-one constraint. Then, based on the obtained optimization variables, this algorithm randomly generates a hybrid beamforming matrix pool satisfying the rank-one constraint by Gaussian randomization method or eigenvalue decomposition method, and selects the local optimum as the approximate solution. However, the performance of the approximate solution of SDR algorithm may be far worse than that of the optimal solution, especially under the highdimensional matrix. To handle this challenge, we invoke the penalty function algorithm [46] [47]. Motivated . .

2) ALGORITHM CONVERGENCE ANALYSIS
This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication. To ensure the effectiveness of the penalty function algorithm, we analyze the convergence of the algorithm. Let the optimal solution of problem 8 Q after the th  iteration be  

F. JOINT DESIGN ALGORITHM OF DIGITAL BEAMFORMING MATRIX AND ANALOG BEAMFORMING MATRIX
In this paper, the final optimization solutions we need to obtain are the digital beamforming matrix and the analog beamforming matrix in the hybrid beamformer, we have obtained the optimization variables   To obtain opt  B , we adopt the eigenvalue decomposition (EVD) algorithm [37], which can be formulated as: where the optimal solution k  b can be given by multiplying the maximum eigenvector of , k opt W by the root of the maximum eigenvalue, the maximum eigenvector and maximum eigenvalue can be obtained by EVD algorithm. Then, we can obtain the hybrid beamforming vectors of K users, i.e., . We can observe that this problem is a joint optimization problem of two matrix variables, which can be regarded as a matrix decomposition problem with power and constant modulus constraints. To handle this problem, we invoke the low complexity algorithms, i.e., NFDM algorithm and AltOpt algorithm [38], which can be formulated as: Note that the columns of the unconstrained optimal digital beamforming matrix are mutually orthogonal in order to mitigate the interference between the multiplexed streams [39]. Inspired by this conclusion, we impose a similar constraint that the columns of the digital beamforming matrix BB W should be mutually orthogonal, i.e., In conclusion, the problem 1 P can be recast as: This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication.

A. COMPLEXITY ANALYSIS
The major complexity of the proposed algorithm is composed of nested iterative optimization and matrix decomposition, and both involve iterative optimization. The number of iterations depends on the predetermined convergence threshold, which is usually small. The complexity of the inner iteration algorithm mainly comes from the CCCP algorithm and the penalty function algorithm. In inner iterative algorithm, problem 8 Q is solved iteratively until it converges to the locally optimal solution. Then, the result obtained by the inner iteration algorithm is substituted into the outer iterative algorithm, to calculate the value of energy efficiency. Next, we update the variable y and continue to next iteration until the outer iterative algorithm converges to obtain the solution. To sum up, the complexity of proposed nested iterative optimization algorithm is approximately I are related to the predetermined convergence thresholds. In the joint design algorithm of digital beamforming matrix and analog beamforming matrix, the complexity mainly comes from alternating optimization in matrix decomposition. In the hybrid beamforming system, the dimension of the analog beamforming matrix is much higher than that of the digital beamforming matrix. Therefor, the complexity of the algorithm is predominated by the analog beamforming part. In each iteration of the AltOpt algorithm, the update of the analog beamforming matrix is realized by a phase extraction operation of the matrix Therefor, the complexity of the matrix decomposition is approximately ( can see that the level of complexity is not high.

B. CONVERGENCE ANALYSIS
In this section, we analyze the convergence of the proposed algorithm. In this paper, we mainly analyze the convergence of the inner and outer nested iterative algorithm, which affects the convergence of the proposed algorithm. In the inner iterative algorithm, problem 8 Q is a typical convex optimization problem, which can be iteratively solved. And the CCCP algorithm is a monotonically decreasing global optimization method, due to the characteristics of the CCCP algorithm, the inner iterative algorithm can be proved to be convergent [43]. For outer iterative algorithm, the optimization problem satisfies the convergence conditions mentioned in Section V-B of literature [33]. In addition, for the AltOpt algorithm, the optimization problem we proposed satisfies the convergence conditions mentioned in [38]. To sum up, the convergence of the proposed algorithm can be guaranteed.

VII. SIMULATION RESULTS
In this section, we illustrate the performance of our proposed robust algorithm in the Massive MIMO LEO satellite communication system through numerical simulations. It is assumed that the variances of the channel errors are identical for different users, which are expressed as 2 Without loss of generality, the SINR constraints are set to randomly generated between 0 5dB : for all users, e.q, . In addition, the values of the adopted the other system parameters are listed in Table II [45]for clarity. This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/ACCESS.2022.3181178  . Meanwhile, we compare the proposed robust algorithm with the traditional non-robust algorithm. In addition, we take the EE under the perfect CSI as a comparison and reference. We can observe from the numerical results that the proposed algorithm can converge to stationary values within very few numbers of iterations, and its performance is better than the traditional non-robust algorithm. And the performance gains become more significant with larger variance of the channel errors.   4 shows the SE performance of the proposed robust algorithm, the channel error parameters are set as above. We can observe from the numerical results that when the channel errors are small, the SE of the proposed algorithm is close to that under the perfect CSI. And performance of the proposed algorithm is better than the traditional non-robust algorithm.   5 shows the change trajectory of EE versus the transmit power threshold T P , the channel error parameters are set as above. We can observe from the numerical results that the proposed robust algorithm is better than the conventional non-robust algorithm under different transmit power threshold. We can see that the EE values initially increase and then decreased as T P increases. This is because the increment of power consumption is faster than that of the system rate.  Fig. 6 shows the EE performance comparison of the proposed penalty function algorithm with Gaussian randomization algorithm and eigenvalue decomposition algorithm [40]. We can observe from the numerical results that the proposed algorithm is better than the other two algorithms. This is because the penalty function algorithm can get a better solution through iterative calculation, which can satisfy the rank-one constraint. Whereas the Gaussian randomization algorithm or eigenvalue decomposition algorithm calculates based on the variables that may not satisfy the rank-one constraint, which may obtain a solution with poor performance.  [41] algorithm and orthogonal matching pursuit (OMP) [42] algorithm. We can observe from the numerical results that the performance of the proposed algorithm is close to the BFGS algorithm and better than the OMP algorithm. The performance of BFGS algorithm is close to the best, which is used as a reference here, but its algorithm complexity is high, and the complexity of each iteration is     . In comparison, although the performance of the proposed AltOpt algorithm is slightly lower than that of BFGS algorithm, its algorithm complexity of each iteration is lower, i.e.,   O MK . Thus, the proposed algorithm has high cost-performance ratio.

VIII. CONCLUSION
In this paper, we have investigated a robust energy-efficient hybrid beamforming design for the Massive MIMO LEO satellite communication system. Taking the CSI errors caused by residual propagation delay and residual Doppler shift into account, we focused on the robust hybrid beamforming design to maximize the system EE under the constraints of transmit power and QoS. Firstly, we adopted the approximate ergodic user rate and the approximate ergodic SINR. Then, we invoked the SDP algorithm to transform the objective function equivalently. Meanwhile, we proposed an inner and outer nested iterative algorithm combining QTFP and CCCP to handle the nonconvex QCQP problem, and adopted the penalty function algorithm to handle the rank-one constraint problem. Finally, we adopted the low complexity NFDM algorithm and AltOpt algorithm to obtain the digital beamforming matrix and the analog beamforming matrix. Numerical results have indicated that the performance of the robust algorithm performs better than that of the conventional one.