The Next Generation of Beam Hopping Satellite Systems: Dynamic Beam Illumination With Selective Precoding

Beam Hopping (BH) is a popular technique considered for next-generation multi-beam satellite communication system which allows a satellite focusing its resources on where they are needed by selectively illuminating beams. While beam illumination plan can be adjusted according to its needs, the main limitation of convectional BH is the adjacent beam avoidance requirement needed to maintain acceptable levels of interference. With the recent maturity of precoding technique, a natural way forward is to consider a dynamic beam illumination scheme with selective precoding, where large areas with high-demand can be covered by multiple active precoded beams. In this paper, we mathematically model such beam illumination design problem employing an interference-based penalty function whose goal is to avoid precoding whenever possible subject to beam demand satisfaction constraints. The problem can be written as a binary quadratic programming (BQP). Next, two convexification frameworks are considered namely: (i) A Semi-Definition Programming (SDP) approach particularly targeting BQP type of problems, and (ii) Multiplier Penalty and Majorization-Minimization (MPMM) based method which guarantees to converge to a local optimum. Finally, a greedy algorithm is proposed to alleviate complexity with minimal impact on the final performance. Supporting results based on numerical simulations show that the proposed schemes outperform the relevant benchmarks in terms of demand matching performance while minimizing the use of precoding.


I. INTRODUCTION
T HE roll-out of the next generation of wireless communication systems is expected to deliver faster internet access and increased capacity, providing customized services in a variety of use cases [2]. Despite the global growth of digital technologies, the United Nations (UN) has recently announced in the General Assembly that half of the world's population still have no internet access [3]. It is because of that there are still many remote locations where fiber and general terrestrial infrastructure cannot be deployed (or not worth the investment), or where the ground equipment is with high probability subject to disruption by man-made events.
Exploiting satellite communications has been identified as a key solution to deliver ubiquitous and high-quality connectivity anywhere in the globe [4]. Conventional High-Throughout Satellite (HTS) systems have employed the spot beam technology with which satellite capacity is equally distributed across the multiple beams and contiguous coverage over a specific region can be provided [5]. While HTS with multi-beam architecture has dramatically improved the satellite system throughput, there have been increasing interests in developing fully reconfigurable satellite schemes that can smartly allocate the high capacity "hot-spot" areas [6].
The recent advances on space technology have opened a door to unprecedented flexibility and adaptability to satellite resources. As highlighted by the major satellite industry experts in Europe [7], "the continuous development of new technologies and the huge increase in satellite interest and investment, witnessed in recent time, have indeed pushed the satellite communication potentialities towards higher limits that need now to be explored to support the efficient and sustainable development of new markets and smart services". In the same document [7], spectrum usage and smart resource management are identified as major research challenges that need to be resolved to unleash the potential of the next generation satellite communication system.
Concerning the satellite industry interest in the aforementioned challenges, next we provide an overview of two of the most advanced GEO HTS systems developed recently. in the world [8]. Coverage, spectrum and capacity can all be reconfigured in-orbit via its innovative reconfigurable payload, to efficiently serve any applications and ensure optimal use of its resources. According to the technical capabilities of Eutelsat Quantum [9], beams can be hopped to spatially diverse regions rapidly and seamlessly. With a similar vision in mind, SES, the satellite operator, worked with Thales Alenia Space to manufacture SES-17. The satellite was launched in October 2021 and incorporates a digital transparent processor (DTP), enabling unique features, such as re-configurable resource allocation, to meet real-time traffic demands [10]. In addition, most of the industry-led projects are still in testing phase, where the algorithm to optimally unleash the flexibility of satellite is still in early stages. For instance, the European Commission has recently launched two 3-year projects related to optimal on-board resource management [11], [12].
Furthermore, similar to the situation in terrestrial domain, the rapid development of data-hungry services has also resulted in spectrum scarcity context in satellite domain. As a consequence, the satellite digital broadcasting standard (DVB-S2X) introduced the possibility to use precoding techniques to enable efficient spectrum management while increasing the spectrum reuse across satellite spot-beams [13], [14]. The feasibility and potential of precoding applied to HTS systems have been recently validated via live experiments on a GEO satellite scenario [15], [16], confirming its relevance for future HTS deployments. In this work, we therefore address a combination of the two aforementioned challenges, namely (i) optimization of payload flexibility; and (ii) spectrum reuse.
Within the flexible satellite payload architectures, this paper focuses on the so-called time-domain flexibility which is commonly implemented via beam-hopping (BH) over time. BH became promising in the early 2010s since this technique can provide a good compromise between complexity and cost. The most attractive feature of BH is the payload mass reduction which is reflected into a reduced cost. Essentially, a BH-enabled satellite scheme can activate a sub-set of beams at each time slot following a time-space transmission pattern and this mechanism can be repeated periodically. In such schemes, the bandwidth can be re-used fully across all activated beams and the interbeam interference can be well-managed by avoiding the geographically-adjacent-beam activation. While BH provides certain degree of flexibility, an extremely asymmetrical traffic demand scenario over the coverage may critically challenge the conventional BH methods. In particular, high-demand areas expanding over multiple adjacent beams necessitate of clusters of beams to be simultaneously activated while making use of the full available spectrum. An example can be the surroundings of an international airport with multiple high-demand mobile platforms flying around or a highly dense populated area with multiple backhauling satellite terminals.

A. Related Works
The works presented in the literature related to multi-beam HTS systems involving BH can be classified into two main categories: 1) Conventional Beam Hopping: The benefits of BH applied to Geostationary (GEO) satellite systems have been well-demonstrated in multiple academic works. Additionally, BH has attracted much attention from some key industrial players, e.g. [17], [18], and has been taken into account in the updated DVB-S2X standard [19]. However, exploiting the conventional BH techniques has also raised some challenges. Conventional BH was conceived to exploit the full available spectrum (i.e. full frequency reuse) over a subset of selected beams, ensuring that the geographical distance between selected beams is far enough to work under a noiselimited scenario [20]. The main technical challenge addressed in the literature has been the design of effective illumination patterns, i.e. determining the different set of beams that need to be activated at each time slot while trying to align the offered capacity with the beam traffic demands over time. The design of illumination patterns for conventional BH usually involves binary variables representing the beam activation simultaneously. Therefore, the problem typically falls within the general mixed integer non-linear programming problem (MINLP) [21], which is very difficult to solve. The authors in [22] considered genetic algorithm targeting a global optimal solution at the expenses of high computational time. In a similar vein, [23] proposed to employ a simulated annealing method which also requires a long time of implementation. As an alternative to optimization-based methods, the works in [24] and [25] developed heuristic iterative procedures which operate in a much faster and more efficient fashion by sacrificing optimality. Clearly, the key challenge identified in early works is the fact that exploiting the beam-activation binary variables results into a large searching space which exponentially aggravates due to the increasing number of potential beams. Following the trends of Machine Learning (ML) applied to optimization problems, [26], [27] investigated the applicability of deep learning tools within the BH illumination pattern optimization procedure. In addition, the conventional BH methods normally focus on no-multiplexing transmission across activated beams, which limits its capability coping with some irregular traffic-demand scenarios in the IoT era [28].
2) Cluster Hopping: The activation of an adjacent set of beams (referred as cluster) was investigated within the European Space Agency (ESA) [28] and proposed in [29] and [30] with the so-called Cluster Hopping (CH) scheme, where linear precoding [14] was considered to mitigate the resulting intra-cluster interference. The downside of the works in [29] and [30] is the fact that the overall virtual multibeam grid is split into a set of non-overlapping clusters of fixed size and shape. This is done to reduce the search space and exclude the cluster design from the optimization problem. While some practical guidelines about the clustering design have been discussed in [30], the problem remains largely unsolved, especially when considering the complexity added by the precoding within the clusters.

B. Our Contribution
In this paper, we propose a general framework for the illumination pattern design, where the transmission of activated beams in separating clusters can be jointly precoded. The objective is to minimize the interference-based penalty with the aim of reducing the use of precoding while constraining the system to satisfy a certain beam demand in a given timewindow. The such technical design can be stated into a Binary Quadratic Programming. Whenever high-demand expands over multiple adjacent beams, the solution from the proposed framework considers precoding to deal with the resulting interbeam interference. With such selective precoding mechanism, complexity at the ground-segment is reduced where precoding operation can be considered flexibly.
To linearize the BQP problem, we first present a procedure to convert the beam demand constraint into a more tractable notation involving required illuminated time-slots per-beam. Next, we present different ways to convexify the BQP problem. First, inspired by the mathematical works in [31] and [32], we reformulate the BQP into a Semi-Definite Programming (SDP) form which can be solved efficiently by employing some standard optimization solver tools. In another approach, we also propose a novel solving framework MPMM which integrates Multiplier Penalty (MP) [33] and Majorization-Minimization (MM) [34] methods. In particular, we relax the binary constraint and add its augmented Lagrangian function into the objective function by using the so-called penalty parameters. Then, the new penaltyform problem is solved iteratively by sequentially updating the penalty multipliers and driving the solution to binary values. In particular, in each iteration, we adopt the MM method to transfer the penalty-form problem into a sequence of simple problems, each of which can be solved optimally. The sequence generated by the optimal solutions of these simple problems is proved to converge to a stationary point. According to the convexity of the penalty-form problem, one also concludes that the stationary point is the optimal solution. Since the previous proposed methods prioritize performance versus computational complexity, we complement this paper by proposing an heuristic greedy algorithm which provides a sub-optimal but efficient solution.
Our main contributions are summarized as follows.
• We propose a general framework and its mathematical formulation to support dynamic and flexible cluster hopping, where geographically adjacent beams are allowed to be simultaneously activated whenever needed according to its demand request. The resulting intra-cluster interference is mitigated with linear precoding, whose utilization is minimized by focusing the design into an interferencebased penalty function. • Based on probability theory, we propose an effective way to reformulate the beam demand constraint and convert it into the number of illuminated time-slots required perbeam in order to satisfy the demand. Such simplification convexifies the constraint and helps easing the tractability of the problem. • Three different methodologies are proposed to address the BQP problem. We first make use of a novel SDP notation specifically designed for BQP problems. As a more accurate alternative, we propose an algorithm that integrates MP and MM methods. Finally, a novel heuristic algorithm is presented to rapidly provide a solution with acceptable performance. • We provide a detailed complexity analysis for each of the proposed methods. • Finally, an extensive numerical evaluation is carried out, where the proposed methods are compared with conventional BH and the previously proposed CH. 1 The results evidence the effectiveness of the proposed algorithms and demonstrate the flexibility of the proposed framework in adapting to any demand distribution. Please note that, although this paper focuses its notation and simulations on GEO satellite systems, the methodology itself can be applicable to the beam placement problem encountered in NGEO constellations. However, the precoding application to distributed satellite swarms is still in early stages of investigation and may need further discussion.
The remainder of this paper is organized as follows. In Section II, we present the system model. In Section III, we present the general formulation of the dynamic beam illumination design problem. To solve the problem, Section IV provide the method to simplify the non-convex demand constraints into the linear forms based on which the BH-design problem is reformulated as a BQP. In Section V, two efficient optimization-based algorithms and a greedy mechanism are presented to deal with BQP. In Section VI, we present numerical simulations. Finally, Section VII concludes the paper. Notations used in this paper are summarized in Table X. II. SYSTEM MODEL This paper studies the forward link of a bent-pipe multibeam geostationary (GSO) satellite system, whose coverage area is divided into a virtual grid of N spot beams. In this system, the illumination pattern is designed over a specific BH window, which is periodically repeated over time. The BH window is divided into a set of M time-slots (TSs), and within each TS no more than K beams (K N ) can be illuminated. The duration of one TS is denoted as T s (seconds) which also represents the minimum dwelling time of the hopping mechanism. Let g l be the traffic demand in bps of beam l, and g = [g 1 , . . . , g N ] T represents the all-beam demand vector. For simplicity, one-user-per-beam scenario is assumed, i.e. a single virtual user located inside the beam footprint (e.g., 4 dB contour) is assumed whenever this beam is activated. Note that the single virtual user is assumed to aggregate the demand of the whole beam user density. The assumption of a single virtual user per beam is performed to abstract the user scheduling. User scheduling is out of the scope of the general BH design for different reasons. The multiplexing of multiple users is assumed to be performed in a time-division-multiple access (TDMA) fashion.

A. Channel Model
Let H ∈ C N ×N be the channel matrix containing all the channel coefficients of the forward link. In particular, the channel between antenna of the satellite payload corresponding to beam l and user k on the ground is modeled following the approach in [5], and can be written as, R is the receiving antenna gain at user k; G l (x k , y k ) stands for the beam pattern gain due to beam l at k which can be estimated according to the user's longitude x k and latitude y k ; φ k,l is the phase component associated with the antenna beam pattern; d k represents the distance from the satellite to that user; λ denotes the wavelength of the carrier frequency band.
Doppler and absorption loss are intentionally not included in our model. The movement of GEO satellite is maintained in a very tight box and has a negligible Doppler shift (note that daily maneuvering is performed to maintain the satellite in its position). Concerning the absorption loss, this would appear as a constant loss in our link budgeting thus not making an impact in our study. This is because the beam hopping window time is usually below few hundreds of frames. For instance, considering a frame duration of 1.3 msec (the number of symbols in a super-frame is 612540, and its duration is about 1.3 msec for a 500 MHz bandwidth), and considering 256 frames, the total BH window-time is approximate to 330 msec. Clearly, the atmospheric loss has a longer coherence time.

B. Selective Precoding Strategy
As depicted in Fig. 1, we envision an illumination pattern design that dynamically activates clusters of beams. Whenever the cluster size is greater than one, precoding is needed to alleviate the inter-beam interference. The illumination pattern design is discussed in Section III. Herein, we detail how precoding is implemented for a given illumination pattern. 2 It is worth noting that the precoding operation entails significant complexity at the gateway side, which exponentially scales with the number of involved beams [35]. Hence, the precoding strategy should be designed smartly by grouping active beams into different precoded-clusters, which are subsequently 2 Note that precoding refers to the exploitation of instantaneous CSI and it is implemented at the ground segment. It needs to be distinguished from the beam pattern formation, which is implemented on-board the satellite and is assumed to be fixed in this work. precoded independently. In what follows, we provide a detail description of the two-step procedure grouping beams in clusters and how the precoding matrix is designed for a specific cluster.
1) Precoded-Clustering Strategy: Given the illumination pattern, at TS t, the definition of precoded-clusters is given as follows. The general idea is that beams generating strong interference to each other should be grouped into one precodedcluster. First, we introduce the so-called influence factor ω i,j , which captures the impact of the inter-beam relative interference from beam i to beam j, and it is defined as where P i denotes the transmission power of beam i and S j stands for coverage area of beam j which is defined by the beam contour at −4 dB from the maximum gain. As can be seen, ω i,j represents the ratio of the interference power from beam i to beam j. Here, the power of signals are calculated as the average value over coverage area of beam j, ( * ) in (2) implies the simplified influence factor when the same transmission power for all activated beams and the same receiving gain at all users are assumed. In a second step, these factors ω i,j are compared to a predetermined threshold. 3 Two beams corresponding to an influence factor greater than the threshold will be located in one cluster. The proposed threshold-based clustering strategy is summarized in Algorithm 1. Particularly, in each TS, one starts by setting every activated beam as a separate cluster. Then, if there are any two beams in two separate clusters that their corresponding influence factor is greater than the threshold, these two corresponding clusters are merged into one. This process is iteratively repeated until there is no change in the clustering structure. Note that the outcome of Algorithm 1 classifies all actives beams into clusters, some of which may contain one beam. Only those clusters with size greater than one will be considered for precoding process.
2) Precoding Design: Once the clusters are formed, the precoding matrix is designed separately for each of them. Let L[t] be the number of clusters in TS t; here, all nonilluminated beams are grouped into a non-transmission cluster for convenience. Denote W i ∈ C ci×ci as the precoding matrix of cluster i where c i stands for its cardinality. For the non-illuminated beams, the corresponding precoding matrix must be zeros since they are silent during that particular TS. For the clusters consisting of only one beam, the precoding matrix can be defined simply as √ P b , where P b denotes the per-beam transmit power. For the remaining clusters (the ones which are formed with more than one active beams), their corresponding precoding matrices are obtained using the MMSE-based strategy as given in [36] and [14]. In particular, the MMSE-based precoding matrixŴ i can be expressed aŝ for l ∈ φ[t] that l > k do 8: if ∃i ∈ C l that max j∈C k ω i,j ≥ κ then 9: where H i represents the channel matrix of all the users in cluster i and α stands for a predefined regularization factor. Then, W i is determined based onŴ i by normalizing every column vector of the matrix to meet the per-beam power ] T including all users' received signal can be expressed as where BDiag( * ) stands for the block-diagonal matrix operation, W[t] ∈ C N ×N represents to the precoding matrix for all beams, s[t] ∈ C N ×1 denotes the transmitted symbol vector; and n stands for the noise vector. In this paper, the zeromean additive Gaussian noise is assumed at all the users where E nn H = σ 2 T I and σ T = √ τT Rx B; τ denotes the Boltzmann constant and T Rx is the clear sky noise temperature of the receiver [5].

III. BH ILLUMINATION PATTERN DESIGN FORMULATION
Let us denote x n,t ∈ {0, 1} as the binary assignment variable indicating the illumination of beam n in TS t. Then, the number of illuminated beams in TS t can be described as N n=1 x n,t . Due to the typical payload mass limitations of a BH-enabled satellite, the number of active beams must remain not greater than the number of RF chains K, i.e., N n=1 x n,t ≤ K, ∀t.
The achievable rate of user n in TS t can be expressed as where w t n denotes the precoding vector designed for user n in TS t, e.g. column of W[t] corresponding to that user. Note that W[t] is also a function of illumination pattern x which determines the way to do clustering for precoding. Considering the average traffic demand of user n, g n [bps], we can express the per-user demand constraint as, The main objective of this work is to develop a BH illumination pattern design such that the users' demands are satisfied while avoiding the use of precoding whenever possible (to ease the complexity burden). For this purpose, we shall avoid the strong cross interference among the illuminated beam in every TS as much as possible. To formulate such problem, we make use of the influence factors defined in (2) to generate a penalty matrix Ω ∈ R N ×N as where [Ω] (i,j) indicates the element on the i-th row and j-th column of Ω. Using this penalty matrix, one can state the BH design problem as 4 where is an integer programming (IP) which is a NP problem in general. The challenge of solving this problem not only comes from the binary assignment variables but also from the non-convex function of R n [t] in constraint (C 2 ) which is a function of x. Remark 1: It is worth mentioning that problem (P 0 ) is stated in a general form in which Ω can take any values. In particular, for different designing goals, Ω can be determined carefully. Therefore, the BH strategy proposed in the following parts can stand in many schemes with appropriate penalty matrices.

IV. PROBLEM REFORMULATION
In this section, effective approaches for dealing with the challenging problem given by (8) are presented. Particularly, according to the idea on estimating the average supplied capacity, we first simplify the complicated non-convex traffic-demand constraints (C 2 ) to linear forms. Based on which, the problem is re-formulated as a binary quadratic programming (BQP).

A. Demand Constraint Simplification
The actual supplied capacity R n [t] is a non-convex function of the variable x t and, therefore, it causes a big challenge for the BH protocol design. Some works in literature have suggested different frameworks to address this issue, such as simple interference-free relaxation given in [37], limiting the set of illuminated beams to avoid the strong interference in [29]. Considering a different approach dealing with this issue, we aim to convert the demand g n [bps] into minimum number of TSs that each beam must be activated in order to meet its traffic demand. To do so, we first consider the following proposition.
Proposition 1: Let ζ n be the average achievable rate of beam n, i.e., ζ n M t=1 . Then, constraint (C 2 ) can be re-formulated as where e n denotes a vector in which the n-th component equals to one and all others are zeros, and * stands for the ceiling operator.
Proof: The proposition can be proved as follows. Since, Additionally, M t=1 e T n x t is an integer. Hence, the right hand side of (11) can be replaced by d n = M gn ζn . This has closed the proof.
In constraint (Ĉ 4 ), d n can be considered as the minimum number of TSs in each of which beam n is illuminated. Next, the following theorem regards the relation between the optimal solution of problem (P 0 ) and d n .
Theorem 1: Let x t 's be the optimal solution of (P 0 ), then we have M t=1 e T n x t = d n , ∀n, if ζ n 's are estimated accurately, which means constraints (Ĉ 4 ) hold for all beams.
Proof: This theorem can be proved easily by using the contradiction method. In particular, one assumes that there exists at least one beam that the corresponding constraint (Ĉ 4 ) does not hold. Denote this such beam as n * which yields M t=1 x n * ,t > d n * . Selecting any TS t * that x n * ,t * = 1, we generate the new solution of (P 0 ), x t 's, that x n,t = x n,t ∀(n, t) = (n * , t * ) and x n * ,t * = 0. It is easy to observe that x t 's satisfies constraints (Ĉ 4 ). Moreover, this new solution meets the requirement of constraints (C 1 ) and (C 3 ) while results in the lower objective function. It follows by a contradict since x t 's is the optimal solution. Therefore, constraint (Ĉ 4 ) holds for all beams with any optimal solution of problem (P 0 ).
Thanks to Theorem 1, the following lemma can stand. Lemma 1: If ζ n 's are estimated accurately, one can replace constraint (C 2 ) by the following one without changing the optimal solution of problem (P 0 ).
Due to this result, in what follows, we propose an iterative framework to estimate the average achievable rate ζ n for each beam by appraising the expected interference.
1) Average Achievable Rate Estimation Framework: According to the mean field theory [38], we assume the uniform distribution of beams to be activated over the time window. Thanks to Lemma 1, the illuminating probability of beam n in a specific TS, e.g., TS t, can be given by Regarding the clustering and precoding processes, the expected interference to beam n can be expressed by taking into account the interference from the beams with low corresponding influence factors and their illuminating probability Based on that, the expected average achievable rate of beam n can be described as Although one of (d n , p n , ψ n , ζ n )'s can be defined if the others are given, determining the accurate values of these factors is very challenging. Exploiting the expectation maximization (EM) algorithm given in [38], we propose an iterative framework to estimate ζ n 's as summarized in Algorithm 2. Particularly, the algorithm initializes with zeros illuminating probability for all beams and repeatedly updating (d n , p n , ψ n , ζ n )'s in each iteration. Each iteration processes two steps, namely, expectation (E-Step) and maximization (M-Step). The E-Step is called for updating the interference based on the illuminating probabilities of the previous iteration while the M-Step stands for calculating ζ n 's, d n 's, and adjusting the probabilities. The iterative process stops at the convergence according to the following proposition. Proposition 2: Algorithm 2 converges after a finite number of iterations.
Proof: As can be observed, p n 's increase while ζ n 's decrease in every iteration. Since, p n 's are upper bounded by ones and ζ n 's are lower bounded by zeros. The iterative process must converge to a stable point after a finite number of iterations.
Remark 2: Note that if the required demand g n is higher than ζ new n in a specific iteration, d new n will be re-set as M -the highest number of TSs. One also notices that problem (P 0 ) is infeasible if K is smaller than the average number of active beams K avg which is expressed as K avg = ∀n d n /M . Remark 3: It is worth noting that (C 2 ) and (C 4 ) may be not equivalent if ζ n is not estimated accurately. In addition, once ζ n is well evaluated as in Algorithm 2, and d n is calculated as in (10), the unmet capacity of beam n must be smaller than, 2) Problem Reformulation: For the sake of simplicity, we compact our notation by rearranging all TSs t into a single tall vector . Thanks to Lemma 1, problem (P 0 ) can be re-stated as Herein, ⊗ denotes the Kronecker product, 1 M stands for the vector with M one elements and I M is the identity matrix with dimension of M .

B. Objective Function Convexification
As can be observed, problem (P 1 ) is a BQP which is NP-hard in general. To ease the tractability of (P 1 ), we aim to characterize the objective function convexity by considering the following theorem.
Theorem 2: For any value of a, problem (P 1 ) is equivalent to the following problem (P a ) : min Proof: Due to the constraint (C 6 ), if x is a solution of problem (P 1 ), then we have ax T Ix = a N n=1 d n which is a constant. Then, x must be a solution of problem (P a ). Inversely, it is easy to prove that any solution of (P a ) must a solution of (P 1 ). Hence, (P 1 ) and (P a ) are equivalent for any value of a.
In addition, the convexity of the objective function (P a ) can be guaranteed if a is selected so that it is not less than −λ(A) -the minimum eigenvalue of A. Thanks to Theorem 2, we can state that problem (P 1 ) is equivalent to a integer QP with a convex objective function, i.e., x TÃ x whereÃ = A−λ(A)I. To this end, instead of solving (P 1 ), we will focus on the following V. BINARY QUADRATIC PROGRAMMING OPTIMIZATION In this section, three optimization approaches are introduced to deal with the BQP problem (P 1 ). Particularly, two efficient solving approaches using the SDP relaxation and MPMM method, respectively. For completeness, a low-complexity greedy algorithm is also proposed. Finally, a complexity analysis for the proposed solution mechanisms is presented.

A. SDP-Based Algorithm
Problem (P 1 ) corresponds to a BQP form, i.e. a problem involving a quadratic objective function with binary variables, which could be solved by relaxing the binary constraint [39]. Firstly, the binary constraint (C 7 ) is equivalent to two equations [40], i.e.
x ∈ {0, 1} MN ⇐⇒ X = xx T and (C 8 ) : diag (X) = x. (22) Herein, the "rank-one" constraint X = xx T can be further relaxed as [40] ⎧ Then, problem (P 1 ) can be approximated to the following semidefinite problem, Problem (P SDP ) in (24) can be solved efficiently by employing the advanced mixed-integer optimization toolboxes such as CVX [41]. If the matrix obtained by solving (P SDP ) is "rank-one", then it provides the optimal solution to the problem. In case of SDP providing a solution matrix whose rank is higher than one, the SDP-based branch and bound method [42], [43] can be applied to obtain the final solution.

B. MPMM Algorithm
In this section, we first introduce the general principles of MP and MM methods based on which a novel multiplier penalty and majorization-minimization (MPMM) algorithm is proposed to solve problem (P 1 ) efficiently. Then, the convergence of this approach is also discussed.

1) Multiplier Penalty Method:
The MP method is an efficient approach for solving the constrained optimization problem. Considering, a general equality-constraint problem as follows, where X is a convex set. Following the MP method, this problem could be solved by minimizing the following sequential problems Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.
where is the index of iteration, {η [] }, {ρ [] } stand for sequences of penalty factors. Here, the penalty term is introduced by its augmented Lagrangian function. The feature of MP method is the way to update η [] step by step [44], which is given by The result given in [33] also concludes that the sequence of {η [] } will converge to a fixed point at which {x [] } converges to a local optimum of problem (25).

2) Majorization-Minimization Method:
The MM method is a well-known approach dealing with a complicated problem by transferring it into a sequence of simple problems which can be solved effectively. The main idea of this scheme is to construct the surrogate function u x | x (k) which approximates the original objective function, then solve the constructed problems in sequence until convergence. For the general minimization problem, min x∈X f (x) where X is convex set. The constructed surrogate function u x Then the sequence of {x (k) } is given by , which will converge to a stationary point of the original problem [45]. If the problem is convex, then the stationary point is the global minimum.
3) Proposed MPMM Method: The challenge on solving P 1 mainly comes from the binary constraint (C7). To cope with this challenge, we aim to employ MP method to relax the binary constraint and deal with a sequence of penalty problems. In particular, the augmented Lagrangian function is added to the objective function with penalty parameters while the binary constraint (C7) is relaxed to form the penalty problem as where f Moreover, the MP-based framework focus on solving (P MPMM ) iteratively and updating penalty parameters η [] , ρ [] to drive the solution of (P MPMM ) to a point that i | is closed to zero. Here, the binary constraint is strengthened by adding 0 ≤ x i ≤ 1, i.e., x 2 i − x i ≤ 0, ∀i to the convex set X . According to the MP-based framework given in [44], the penalty parameters can be updated as where β [] is the parameter to update ρ step by step. Usually, ρ [] would be initialized with a small value and then increase with iterations. Note that ρ [] can also keep fixed after certain iterations [33]. Due to the high order of variable x appearing in the objective function, it is very challenging to solve (P MP MM ) directly. To overcome this issue, we employ the Algorithm 3 MPMM Algorithm 1: Initialization: = 1, η [1] = 0, ρ [1] = 1, x [1] (0) = 0 2: repeat 3: Set k = 0. 4 MM method by constructing surrogate function and then finding the optimum solution in a sequence. The proposed algorithm is summarized in Algorithm 3. The surrogate function is constructed by linearizing the binary constraint with its first order Taylor series and is given by

4) Convergence Analysis:
Regarding the same approach analyzing the convergence of the MP method given in [33], one can prove the convergence of the proposed algorithm by addressing two facts: i) for given η [

C. Greedy Algorithm
In this section, we propose an heuristic approach solving (P 1 ). The basic idea is to firstly solve the relaxed problem, i.e. (P rlx ), in which the binary variables are relaxed as continuous one, i.e.
The continuous optimal out-comes, denoted as x con , are then rounded to binary solution, x bin . The rounding mechanism is developed so that all the practical requirement of (P 1 ) are guaranteed. In addition, the approach also aims to minimize the total penalty. In particular, after solving the problem (P rlx ), for each user n, the first d n highest elements among the set {x con n,t } M t=1 are set to be ones while the others are down-rounded x bin j,t , ∀(n, t). 3: while x bin is not a feasible solution do x bin n,t < K . 5: Let T = (n, t, u)|t ∈ S + , u ∈ S − , x bin n,t = 1, x bin n,u = 0 . 6: Solve (n * , t * , u * ) = arg min (n,t,u)∈T B n,u − B n,t . 7: Swap the values of x bin n * ,t * and x bin n * ,u * by setting x bin n * ,t * = 0, x bin n * ,u * = 1. 8: end while to zeros. Then, the illumination of beams over TSs can be further swapped to satisfy constraint (C 5 ) and keep the total penalty as small as possible. The greedy algorithm presented in Algorithm 4 and described in Fig. 2.

D. Complexity Analysis
In this section, the complexity of our proposed algorithms is investigated based on the number of required operations.
1) Threshold-Based Clustering Algorithm: Let B t = |B ac [t]| be the number of activated beam in TS t where |X | stands for the cardinal number of set X . As can be observed, Algorithm 1 consists of three loops, i.e., one "repeat" loop and two "for" loops, and it initializes with B t clusters each of which contains one beams. In each iteration of the "repeat" loop, the number of cluster decreases if Φ[t] changes. Moreover, the "repeat" loop stops when there is no change of Φ[t]. Therefore, the iteration number of "repeat" loop must be less than B t . Regarding "for" loops, we can observe that for each couple (k, l) in φ[t], one has to compare ω i,j to κ for all (i, j) ∈ C l × C k . Then, according to |φ[t]|, |C k | ≤ B t for all k ∈ B ac [t], the complexity of Algorithm 1 can be estimated 2) SDP-Based Algorithm: As given in [46], the computational complexity involved in solving the SDP is O(max(m, n) 4 n 1/2 ) where n and m are the numbers of variables and constraints, respectively. As can be observed, these numbers corresponding to problem (P and I inner MPMM are the average iteration numbers of outer and inner loops required in Algorithm 3 to solve problem (P 1 ), respectively.
4) Greedy Algorithm: The initial step of this algorithm attempts to solve the QP (P 1 ) with continuous variable x with the complexity of O (M N ) 3 [47]. Then, it requires to calculate elements of matrix B before iteratively updating sets S + , S − , and T . At the end of each iteration, the comparison procedure is processed to select the two specific beams in two different TSs for swapping. It is worth noting that the cardinality of T cannot exceed M N and decrease after every iteration. Hence, the number of iterations in Algorithm 4 must be smaller than M N . Hence, the complexity of the greedy algorithm can be estimated based on that due to solving QP, calculating B and iterative process as O

VI. NUMERICAL RESULTS
Monte Carlo simulations are conducted to evaluate the performance of the proposed three algorithms in terms of demand matching, the total computation for precoding and average number of active beams per TS, comparing with two benchmarks: conventional beam hopping (BH) and cluster hopping (CH). In particular, the conventional BH method is given in Appendix A while CH benchmark solution is proposed in [29] and [30].

A. Simulation Setup
We consider a GEO satellite system with 67 spot beams, i.e., N = 67. The setting parameters are summarized in Table I. The simulation setup is the same as the one considered in [28] and [30]. Unless mentioned otherwise, the number of TSs is set to M = 20. The traffic demand of all the users are generated uniformly at random between 400r and 1500r (Mbps), i.e., 400r ≤ D n ≤ 1500 r ∀n. Herein, r represents the demand-density factor which is selected in {0.25, 0.3, 0.35, 0.4, 0.45} where r = 0.25 implies the low  demand setting while r = 0.45 refers to the high demand. For each selected value of r, 50 demand instances are generated for testing. A single representative user within each beam is assumed, which aggregates the overall beam demand. the users' demand cannot be satisfied if K is less than K avg ; hence, unless mentioned otherwise, we set K = K avg in the subsequent simulations.

C. Discussion on MPMM Convergence
1) The Convergence of MPMM Algorithm: In Fig. 4, we regard the convergence of Algorithm 3. In order to illustrate the convergence, three parameters are considered, i) gap (x i ) = min{|x i − 0|, |x i − 1|} which describes the minimal distance between the continuous element x i and a binary variable; ii) T (z) = {x i | gap (x i ) ≤ z, i = 1, . . . , MN } is the set of elements which belongs to the variable x, whose gap (x) is not greater than z; iii) P (z) = |T (z)| MN describes the percentage of the elements in variable x whose gap (x) is less than z. Fig. 4a shows the geometric distribution of elements in x, i.e. P (z) achieved in each outer-loop iteration. For the sake of clarity, in this simulation a single demand factor r is considered, being r = 0.3. It can be observed that the elements in x are closer to binary values after each iteration. As can be observed, the curves corresponding to iteration 3 and 4 illustrates that P (z) close to one with very small gaps. Certainly, this has confirmed that the final solution converges to binary variables. Fig. 4b shows the variation of the objective function of problem (P 1 ), i.e., x TÃ x, achieved in every outer-loop  iteration. The algorithm initiates with penalty parameters η [1] = 0, ρ [1] = 1 which will increase after each iteration.

D. Performance Evaluation
In this section, we aim to evaluate the performance of the proposed algorithms in terms of two main aspects: (i) perbeam demand matching, and (ii) number of beams that would require implementation of precoding to deal with co-channel interference.  Fig. 5a shows an example of an illumination pattern design obtained by implementing Algorithm 3 for a particular demand instance obtained with r = 0.25. In this figure, the white rectangles imply that the corresponding beams are illuminated while the blacks refer to the inactive ones in a specific TS. In addition, illumination map of beams corresponding to TS 15 of this simulation is illustrated in Fig. 5b, where the green areas represents the foot-prints of the illuminated beams.
Next, we consider the precoding utilization in the proposed algorithms in Table II, III and IV, which will have an impact on the system complexity. In is worth noting that the complexity for MMSE-based precoding of a cluster of N beams is estimated as O N 3 in general [14], [36]. Therefore, we aim to demonstrate the precoding complexity corresponding to different BH mechanisms by illustrating the number of clusters with different sizes. In particular, Tables II, III and IV show the distribution of clusters in various sizes at different demand instance assuming r = 0.25, 0.35 and 0.45, respectively. In addition, based on the numbers given in these tables, the total computation for precoding due a specific BH mechanism can be estimated as T = N ni=2 m i ·(n i ) 3 where n i represents the size of clusters and m i represents the total number of the corresponding clusters at a demand density instance. At footnote of these tables, we show the total relative computation cost for precoding by comparing T 's, where SDP is set as the baseline. Particularly, the precoding complexity due to a BH method is defined as the ratio of its T to that of SDP. For the three proposed algorithms, these tables have demonstrate that: (i) the number of larger-size clusters is smaller, (ii) increasing the traffic demand results in the higher number of larger-size clusters. Interestingly, MPMM method shows its superior when it on the smaller size of cluster which would result in less computation for precoding. When r = 0.25, the total computation for precoding with SDP method is more than 300 times of that of MPMM. In addition, no adjacent beams will be illuminated simultaneously with BH method, then there is an upper bound of the maximum number of activated beams. So the total number of active beams will not change much with the increase of the density of demands. The CH method predefines the clusters where the cluster size is 6 or 7.
Next, we analyze the average number of active beams per TS, which determines the resulting interference as well as the operating power consumption of the satellite. In principle, one would like to minimize the number of active beams but making sure that the demand requirements are met. Table V shows the average number of beams activated per TS for the different methodologies. For the proposed algorithms, the number of activated beams per TS is fixed and given by Algorithm 2, while the conventional methods provide different values. For the CH technique, the number of active beams is fixed and does not depend on the demands, which typically results in an inaccurate demand-matching performance. Regarding the BH technique, the number of active beams slightly increases as the demand increases, but the illumination design is limited to non-adjacent beams, and therefore the increase in number of active beams is not so prominent. Unlike the benchmarks, the proposed techniques are more flexible in activating more number of beams and adapting to the demand increases. The results in Table V match the distribution of of cluster number  with different size depicted in Tables II, III and IV. To evaluate the fairness of users' satisfactory corresponding to the proposed and benchmark methods, we consider the Jain's Fairness Index proposed in [48]. The definition of the index is given as J (y) = ( where y i is the chosen metric and is given by y i = ci gi in which c i is defined as . This index aims to determine whether users are receiving a fair demand matching or not. "One" value of J (y) implies the highest fairness level among all users. The Jain's fairness indices achieved by implementing various BH mechanisms, the proposed and benchmark methods, for different demand factors are given in Table VI. From the results, it can be concluded that all the three proposed methods can provide better fair indices than the benchmark, specially in the high-demand scenarios. In addition, it can be observed that the CH method stays at around 0.86 independent of the demand entry. Furthermore, the Jain's index of the conventional BH method suffers to maintain a good level of fairness as the demand increases.
In the subsequent results, we focus our evaluation on the capabilities of the proposed techniques to match the offered capacity with the actual demand. In particular, Fig. 6 illustrates the cumulative distribution function (CDF) of C/D -the ratio of provided capacity of the beam to its required demand, for 3 different demand factors r = 0.25, 0.35 and 0.45. On the top of Fig. 6, we show the performance of the proposed methods with respect to the benchmarks. On the bottom of Fig. 6, the reader can find a zoom-in figures to better discern the performance of the proposed techniques. The vertical dashed line indicates the ideal scenario where c i /g i = 1, ∀i.
First, Fig. 6 confirms that the CH technique suffers from the limitation of pre-defined clustering shapes, which unavoidably illuminate low-demand beams with high-demand beams. On the other hand, the conventional BH is shown to experience significant degradation when the demand factor increases. This is because BH falls short in supplying enough capacity due to its inability to illuminate high-demand areas at once. Focusing on the proposed techniques, we can see from Fig. 6(b) and Fig. 6(c) that all three outperform the benchmarks, specially for moderate and high demand factors. It can also be observed that SDP-based method provides c i /g i > 1 for almost all cases which implies that this approach can provide the capacity larger than the demand. It may not be expected in some specific circumstances which one avoids spending expensive network resources to serve the users much more than what satisfies them. The MPMM approach seems to provide a better trade-off as its curve is closer to the ideal case. The greedy algorithm provides a performance in between SDP and MPMM methods, and seems to be closer to SDP solution for low demand factor while it approaches the MPMM solution for the high demand factor.

E. Impact of Number of TSs in BH Window
Herein we evaluate the impact of the parameter M , which determines the number of TSs within a BH window. In particular, Fig. 7 evaluates the number of average precoded beams within a hopping-window with respect to M shows the CDF of (c i /g i ) for different values of M . In Fig. 7, we focus on the MPMM method's behaviour, which was found to be the best in demand matching fairness among users in the previous simulation result. As usual, we evaluate three different demand factors, i.e. r = 0.25, 0.35 and 0.45. From Fig. 7(a), Fig. 7(b) and Fig. 7(c), it can be observed that a higher value of M translates into a lower average number of precoded beams. This is an expected result where as the more TSs are available, the less number of beams need to be activated simultaneously. Another interesting result is that the average number of precoded beams increases with the users' demand. Focusing now on the CDF curves, depicted in Fig. 7(d), Fig. 7(e), and Fig. 7(f), we can observe that the longer window length can push the achievable capacity  closer to the users' demand. This is because of that the higher value of M gain the higher degree-of-freedom in selecting the illuminating TSs for each beam to meet their demand.

F. Impact of Imperfect CSI
Simulations above based on perfect CSI which is unrealistic in practical. In this subsection, following [49], we model the CSI uncertainty with an additive complex Gaussian error with parameters (mean, standard deviation) as shown in Table VII. Note that (I/N ) in Table VII denotes the Interference-to-Noise Ratio, which is a measurement of how strong are the signals (coming from different beams) to be measured. Table VIII compares the performance with perfect CSI and imperfect CSI in terms of estimated demand (d n ) and Jain's fairness index for different demand densities r, where the results are averaged for 50 Monte Carlo simulation. The only difference between each simulation is the channel information, one of which is with perfect channel and the other is with estimated channel. First thing we observe is that the estimated average demand d n is lower with imperfect CSI. The latter occurs due to the nullification of certain CSI components in the imperfect CSI (note that I/N lower than −10dB are not measured at all). The nullification of the channel matrix translates in a reduction of the assumed interference levels. This biased demand estimation could be compensated by modifying the way we calculate the average demand, maybe adding a margin, but this is out of the scope of this work. When comparing the values of Jain's fairness index in Table VIII, we can observe an evident performance loss for the imperfect CSI case, which is justified essentially by the reduced estimated demand. Fig. 8 illustrates the demand matching for the same cases evaluated in Table VIII for completeness. As expected, the figure shows that the perfect-CSI scheme outperforms the imperfect-CSI one where it can supplies more beams as their demands than the other.

G. Impact of Random User Location
The assumption of a single virtual user per beam is performed to abstract the user scheduling. However, the assumed location of such virtual user may have some impact on the final    performance. For instance, having the virtual user on the beam edge will not have a strong impact whenever the active beam is isolated. However, for high-demand instances, we expect the edge users to impact on the selective precoding and generate a higher miss-match between the estimated capacity and the actual supplied capacity. For the latter, estimating the capacity with a user on the edge will provide lower capacity than a user in the beam center, therefore requiring more number of TSs to satisfy the demand.
To evaluate this, we have run some results by randomly selecting the virtual user location within its −3 dB beamwidth. For the sake of comparison purposes, for each instance of random user location, the same beam traffic demand (in bps) as the user in the beam center is assumed. In addition, there are 50 instances, each of whose demands are randomly generated. Table IX compares the performance in terms of demand satisfaction between random user and centered user at different density of demand. It can be observed that, as the demand r increases, the error in the demand matching increases when non-centered virtual user is considered.

VII. CONCLUSION
In this paper, we propose an analytical framework, a class of BQP problems, to support dynamic beam illumination design considering selective precoding for the next generation of time-flexible satellite broadband systems. Three algorithms are proposed to solve the problem: (i) SDP-based approach, (ii) MPMM methodology and (iii) low-complexity greedy algorithm. All three methods target the cross-beam interference minimization, such that the number of beams that need to be precoded are kept to minimal in an attempt to reduce system complexity.
An extensive evaluation has been carried out based on numerical simulations. The results have shown interesting gains provided by the proposed algorithm with respect to the relevant benchmark schemes. In particular, the proposed framework provides an efficient solution to deal with highdemand areas while keeping the precoding-related complexity low.

APPENDIX A CONVENTIONAL BEAM-HOPPING METHOD
The conventional BH method is one of the methods to design the illumination pattern and is developed by solving the following problem where X = [x 1 , x 2 , · · · , x M ] and B = {(i, j) |Q i,j = 1, i = 1, · · · , N; j = 1, · · · , N}, Q is the adjacent matrix of the graph G = (V, E) and ζ ∈ R N represents the estimated capacity for all the beam, which is given by (14) where no inter-beam interference is considered. Herein, G = (V, E) is defined as follows. Each beam center is considered as a vertex v ∈ V and any two vertices are connected with an edge e ∈ E if those vertices represent geographically adjacent