Inbound Carrier Plan Optimization for Adaptive VSAT Networks

The past decades witnessed the application of adaptive coding and modulation (ACM) in satellite links. However, ACM technologies come at the cost of higher complexity when designing the network's carrier plan and user terminals. Accounting for those issues is even more important when the satellite link uses frequencies in Ka band and above, where the attenuation caused by tropospheric phenomena is a major concern. In this article, we propose a solution for the inbound, i.e., return link, carrier plan sizing of very small aperture terminal networks. As tropospheric attenuation is a key factor, we present a mathematical problem formulation based on spatially correlated attenuation time series. Our proposed sizing scheme is formulated as a mixed integer linear programming optimization problem. The numerical results for a test scenario in Europe show a 10 to 50% bandwidth improvement over traditional sizing methods for outage probabilities lower than 1%.


I. INTRODUCTION
For satellite communication (SatCom) systems, radio bandwidth is one of the most precious resources. Satellite payloads are designed to deliver services for a variety of applications, each associated with a target service level agreement and desired area of service. The amount of bandwidth required to deliver services while meeting the service level agreements can vary considerably, depending on the type of service provided, such as TV broadcast, video on demand, voice, and Internet. Since spectrum is limited by regulations, SatCom service providers strive to deliver services to their customers using as little bandwidth as possible. Bandwidth optimization is, thus, of paramount importance.
When optimizing bandwidth usage, we can distinguish two scales.
1) Micro/real-time scale: In this case, the optimization focuses on real-time bandwidth and time-slot allocation [1], with short-term impairments and traffic prediction [2]. Techniques such as connection admission control [3], channel estimation [4], power allocation [5], and adaptive coding and modulation (ACM) [6] fall into this category. 2) Macro/system scale: In this case, the optimization is done at a system level, dealing with long-term statistics of channel behavior, link budget outputs, and traffic variations. Optimization on this scale usually aims to reduce the cost per gigabit per second [7] by tweaking key satellite characteristics, such as beam number, frequency reuse, and signal-to-interferenceplus-noise ratio (SINR) [8], or using smart gateways and temporal/spatial availability [9].
One of most well-known system-scale optimization problems in the current literature is the throughput capacity optimization of both high throughput satellites and upcoming very high throughput satellites, with flexible payloads (especially for the forward link) [7], [10]. This kind of optimization works best before the satellite deployment when the payload can still be modified ("Upstream" part of Fig. 1 [11], [12]).
The topic of this article also falls into the system-scale category but focuses on solving a different problem when the satellite is already deployed. Once in orbit, the satellite's transponders will be used to provide space segment resources, such as bandwidth, to support service provider's broadband networks, each using a fraction of the satellite's spectrum (service providers are located in the "Downstream" part of Fig. 1). Each service provider's request will give rise to a different network design (or modification of an existing one), which will require its own optimization process. The network designers also have to take into account specific constraints, thus increasing the complexity of the optimization problem. Limitations currently found in the industry most notably include the following.  [13]. In this case, the return carriers are not allowed to adapt their ModCod or bandwidth to the fade conditions-user terminals can only select one of the pre-existing carriers. It is then necessary to establish a carrier plan, best described as a matrix C, with its entries c k,l being the number of carriers associated with M k and B l .
Careful design of the carrier plan is an extremely important step of the design process. On the one hand, the carrier plan dictates the bandwidth B tot used by the network via the following equation: On the other hand, the carrier plan also dictates the throughput capacity of the network R tot via the following equation: where S k denotes the spectral efficiency of M k in b/s/Hz. Therefore, the carrier plan must also be designed to deliver the contractual quality of service (QoS) to the terminals, which can conflict with the bandwidth minimization. The design of the carrier plan is then a complex tradeoff between bandwidth minimization and QoS satisfaction. The complexity is at its highest when designing the return link sizing of very small aperture terminal (VSAT) networks since: 1) the link budget is tighter in the return link compared to the forward link, thus the carrier types (ModCod and bandwidth size) are more heterogeneous in order to offer more options in case of fade; 2) a potentially large number of terminals with stringent QoS requirements may be considered. The primary goal of this article is to determine what the optimal VSAT carrier frequency plan needs to be in order to both meet the QoS constraints and minimize the total bandwidth consumption. To the best of our knowledge, the problem of optimal inbound carrier plan design, and more generally return link sizing, in the adaptive VSAT context has not been satisfactorily addressed in the scientific literature nor in the industry. On the industry side, service providers use heuristics without optimality guarantees [14]- [16].
In the literature, the sizing of fade-adaptive satellite networks has been studied in [17] to demonstrate how time-series synthesizers can improve ACM network designs. Luini et al. [17] have generated a time series of rain maps using the MultiEXCELL rain field simulation model and have used them to determine the ModCod usage statistics. Using a long-term average approach, the authors have presented how to derive the necessary number of carriers and have shown how to build the complementary cumulative distribution functions (CCDF) of the number of carriers needed for each ModCod. However, one of the key questions left unaddressed by the authors is how to use such CCDFs to find the optimal carrier plan.
Pranjic et al. [18] have tried to answer this question. They have focused on using a mixed integer linear programming (MILP) formulation to optimize the carrier frequency plan, allowing to take into account hardware constraints such as a limited set of carrier bandwidths. However, the resulting MILP formulation falls short of providing a practical solution due to its inability to guarantee the QoS. The proposed solution also fails to take into account the spatial correlation of terminal's impairments, as the methodology focuses on individual terminal CCDFs and does not introduce correlation matrices. This work has inspired our previous work on optimizing inbound carrier plans for constant coding and modulation networks [19].
In this article, we expand the works in [17]- [19] with the following new contributions. 1) A formal definition of the return link carrier plan sizing optimization problem using spatially correlated attenuation time series from the input definitions to the identification of optimization variables and constraints. 2) An MILP problem formulation which is solvable by standard solvers, as well as a comparison against classical sizing methods in a test scenario in Europe using MultiEXCELL for the generation of tropospheric attenuation.
The rest of this article is organized as follows. In Section II, we introduce the VSAT ACM scenario as well as the carrier plan definition. In Section III, we define the attenuation time series along with other input parameters and show how to formulate a general optimization problem. In Section IV, we present a reformulation of the optimization problem, solvable by MILP solvers, as well as a In Section V, we introduce a test scenario in Europe, and we present numerical results comparing our MILP method to the benchmark. In Section VI, we discuss some aspects of the two methods and the corresponding results. Finally, Section VII concludes this article.
Notation: Bold uppercase C denotes matrices, bold lowercase x denotes (row) vectors, and calligraphic uppercase M denotes sets. |M| refers to the cardinality of set M, and · refers to the ceiling operator. For the reader's convenience, Table I describes the key notations found in this article.

II. SCENARIO
This section describes the context and the different assumptions made in this article. We first present the carrier plan and the key constraints it is subject to. We then present simulators of spatially correlated tropospheric attenuation and explain why they are relevant to the return link sizing problem.

A. VSAT Inbound Carrier Plan
We consider the return link carrier sizing problem of a single-beam satellite network serving N VSAT terminals. We assume that all terminals require the same uplink services. We also consider that the following QoS inputs have been agreed upon or computed in previous design steps.
1) The committed information rate (CIR) R CI (in kb/s) which is the contractual minimum bitrate under which a terminal is considered in outage.
2) The target outage probability p, which is the maximum fraction of time a terminal can be allocated a throughput lower than R CI . It is mainly dependent on the fade margin of the link and the robustness of the carrier plan to congestion.
Our goal is to find a carrier plan that minimizes the network's required satellite transponder bandwidth while ensuring that the QoS constraints are satisfied. Given a pool of ModCods available in the network M = {M 1 , . . ., M K } and a pool of bandwidths B = {B 1 , . . ., B L }, we defined in Section I a carrier plan as a matrix C, with c k,l being the number of carriers associated with M k and B l . Since, in [19], we have proposed a heuristics that, starting from a given proportion of terminals x k ∈ [0, 1] using ModCod M k , allows to compute the optimal c k,l values, for all k and l, that minimize the total carrier plan bandwidth, then we will focus on computing the optimal x k in this work. We interpret x k as the maximum proportion of terminals the network will be able to serve the QoS with ModCod M k . We then define the (simplified) carrier plan x as As the network uses ACM schemes, terminals will be able to request different ModCods over time. The request of a ModCod by a given terminal will be driven by the attenuation conditions at this terminal's location. Fig. 2(a) illustrates an example of the terminals' requests with N = 4 terminals and K = 3 ModCods. The yellow terminal is in clear sky and can use the most efficient ModCod. On the other hand, the blue and green terminals, experiencing moderate fade, can only close a link using more robust ModCods than that available in clear sky. Finally, the red terminal is facing severe fading and cannot close the link.
The carrier plan, however, will not be allowed to change over short periods of time. This configuration, also used in [18], is inspired from the ARCS technology. In Fig. 2, we also illustrated an example of such carrier plan with . This carrier plan can then serve at best when one terminal using M 1 , one terminal using M 2 , and two terminals using M 3 for a total of four terminals. We represented the terminal capacity for each ModCod with small boxes, where one box represents a capacity of one terminal, and the color inside the box represents the terminal using this capacity at this time instant. For instance, in Fig. 2(a), the blue terminal is using M 1 , the green terminal is using M 2 , and the yellow terminal is using M 3 . Only one terminal is able to use M 3 , therefore half of the capacity provisioned for M 3 is unused. This unused capacity is represented as a blank box. With this approach, our work aims to find practical solutions to common cases in the industry where the assumption of fully flexible carrier plans does not hold because of technology limitations.
In this context, great care must be taken in the composition of the carrier plan to avoid congestion, i.e., to prevent too many user terminals requesting the same ModCod at the same time. It is expected that, given the spatial correlation characteristics of the fade, it is likely that at a given moment in time, certain carriers will experience higher demand than others due to the number of terminals experiencing similar levels of total link degradation and only able to close links using those specific ModCods. Beyond a certain congestion level, terminals will be denied service. Fig. 2(b) illustrates such a case, as the blue and green terminals, facing the same fade conditions, have to request M 1 . However, the carrier plan x, identical to Fig. 2(a), can only afford one terminal using M 1 , and consequently, the blue terminal is dropped due to congestion. Furthermore, the acceptable proportion of time service is denied to a given terminal is limited by the outage probability p, as defined earlier. Thus, the spatial correlation of fade plays a critical role in the design of the carrier plan.

B. Simulators of Spatially Correlated Attenuation
Over the past decades, SatCom systems have explored unoccupied spectrum regions wherein wider bandwidths are available. Many recent satellites use Ku and Ka bands, and more tests are being conducted using the Q/V band [20], which are expected to be employed by future telecommunication satellites [21]. However, an increase in frequency also leads to an increase in the tropospheric attenuation of the radio signal, thus requiring more accurate models to enable a proper design of such communication systems [22].
While most generally accepted models have been integrated into the well-known ITU recommendation propagation series [23], the statistical models have been criticized [24] for not being designed to make predictions with high spatio-temporal resolution. To address this issue, timeseries synthesizers of tropospheric fade, considering both time and spatial domains, have been developed [24]- [26]. While, first, only able to generate attenuation for short periods (few hours) and small scales (hundreds of square kilometer), newer models such as in [27] have increased the size of the space and time dimensions of the synthesized series, allowing continental scale simulations spanning multiple years. These breakthroughs have fostered the use of time series generators in the design of SatCom networks [28].
The application of spatially correlated time-series synthesizers has been first explored for the forward link in [29]. This work has illustrated how time series generators could be used to compute network-scale statistics, key to enabling applications in network sizing. Joint terminal fade statistics have been used to compare ACM and variable coding and modulation (VCM) technologies in a multibeam video broadcast satellite scenario [30]. Their numerical simulations with the software SISTAR have shown that VCM performs better in the context of multibeam forward link broadcast.
The impact of fade correlation was found to be even more important on the inbound side. The impact on the Shannon capacity of multibeam return links has been thoroughly explored in [31]. The authors have used a statistical model of the rain fade correlation, similar to [26], in the modeling of the terminals' channels. From this model, the authors have shown how to deduce long-term CCDF of the network outage capacity in function of the outage probability and have concluded that introducing spatially correlated rain attenuation results in a significant decrease in channel capacity.
In order to consider the spatial correlation of fade, we propose, in this article, to take the time series of spatially correlated attenuations as an input of our design method.

III. PROBLEM FORMULATION
The outputs of a time-series synthesizer are the attenuation time series a n (t ) faced by the terminal n ∈ {1, . . . , N} at time t ∈ {1, . . . , T }. Knowing the satellite and ground stations characteristics, the SINR C N+I n (t ) can then be computed through a link budget as follows: is the carrier to intermodulation noise ratio, and (5) and where 1) EIRP 0 n is the equivalent isotropic radiated power (EIRP) density of terminal n. We assume the terminals work in power equivalent bandwidth (PEB) mode in which emission power is scaled with bandwidth.
2) L fs n is the free space loss between terminal n and the satellite. Along the C N+I n (t ) time series, we assume the following additional parameters are given.
1) The number of terminals N.
2) The CIR R CI in kb/s. The output of our optimization process will be the optimal carrier plan x opt ∈ [0, 1] K , which minimizes the system bandwidth under the constraint that no terminal will be delivered a throughput inferior to R CI for a proportion of time greater than p. For the sake of conciseness, we will refer to this constraint as the "outage constraint."

A. ModCod Request Time Series
Before diving in the optimization problem, we first need to process the SINR time series in order to extract a time series of ModCod requests, which we can then use to find the optimal carrier plan. An illustration of the different steps is given in Fig. 3. After computing the SINR time series from the attenuation, the following step is to find the index K n (t ) of the most spectrally efficient ModCod M K n (t ) that the terminal n can use at time t. This is done by comparing the terminal SINR, C N+I n (t ), against the ModCod thresholds in order to satisfy the following: We can then determine the ModCod request proportion time seriesx k (t ) as follows: for all k ∈ {1, . . . , K}, and In this article, we will use the hat notation to indicate that x k (t ) is a sample computed from the SINR time series. We also define λ n , the proportion of samples in which terminal n is not able to close the link, as We assume that, for any terminal n, λ n is strictly less than p. Having this condition violated would mean the terminals have not been designed properly in early design steps, and no carrier plan will be able to satisfy the outage constraint. From a given x, we are able to compute the bandwidth needed to deliver R CI to each of the terminal proportions x k via Our goal is then to find the optimal carrier plan with O ⊂ [0, 1] K the viable set of carrier plans that satisfy the outage constraint, which we need to express mathematically.

B. Outage Constraint Through Order Relation
Let us assume we have found x opt and have selected it as our network carrier plan. Let us also assume that our network terminals meet fade conditions that lead to a ModCod demand equal tox(t ). In this scenario, the proportion of terminals requesting M k isx k (t ). If x opt k is greater or equal tox k (t ), all the corresponding terminals will have their demand satisfied. However, if x opt k is less thanx k (t ), the network's carrier plan will not be able to satisfy all terminals' demands. The remaining proportion of terminalsx k (t ) − x opt k will have to seek free capacity in other coordinates of x opt , i.e., use another ModCod. As ModCod usage is conditioned by the terminal SINR, the remaining proportion of terminals can only be reallocated to x opt i with i ∈ {1, . . . , k − 1}. Having a capacity allocated to every terminal requesting M k at time t, including the remaining terminals, is then equivalent to having at least x k (t ) − x opt k free capacity in the k − 1 ModCods between M 1 and M k−1 , i.e., By iterating on k, we deduce that all terminals able to close the link at time t can have a capacity allocated in We then denote the order relation as (14) and say that x opt can accommodatex(t ). This relation is a partial order, but it is not a total order, as it is possible to have both x opt x(t ) and x opt x(t ), and thus, x opt x(t ) x opt x(t ). It is also different from the wellknown majorization relation [32] since the vector's coordinates are not sorted in descending order.
Let A x denote the set containing all samples accommodated by x ∈ [0, 1] K as follows: If terminals could always close the link, we could ensure the satisfaction of outage constraint by designing x opt such that no terminals will be dropped because of congestion for at least T · (1 − p) samples. This is equivalent to have x opt accommodating T · (1 − p) samples and is written as Note that (16) is more stringent on x opt than the original outage constraint. It implicitly states that x opt cannot drop even one terminal for at least T · (1 − p) samples. For the remaining T · p samples, one or more terminals will be in outage but it will likely not be the same terminals for every sample. Thus, even the highest outage statistic among terminals can be much lower than p in certain cases. However, the terminals outages are not just due to congestion but also due to link outage. Let λ max be the maximum proportion of time during which one terminal cannot close the link, that is As we do not know when a terminal will have its link in outage, we need to assume that the set of time samples at which a terminal cannot close its link is disjoint from the set of time samples at which this terminal is dropped due to congestion. This implies that x opt will have to accommodate not just T · (1 − p) samples but T · (1 − p + λ max ) .
The optimal carrier plan then has to satisfy We can finally express O as and, therefore, rewrite (11) as the following optimization problem:

IV. OPTIMIZATION FORMULATION
The main difficulty in solving (20) resides in efficiently determining if a given x satisfies (20b). In this section, we reformulate (20b) with integer constraints, enabling the use of MILP solvers. Then, we introduce a worst-case method that we will compare against our MILP formulation in Section V.

A. Mixed Integer Linear Programming
The optimization problem (20) is quite difficult to solve due to the cardinality constraint on A x in (18), which is nonlinear. Drawing inspiration from the work in [33], we can translate the cardinality constraint into an integer constraint. This is done by introducing a binary optimization variable, α ∈ {0, 1} T , as in (21) shown at the botom of this page. This variable lets the solver decide which samples should be accommodated by x. We can then rewrite the cardinality-constrained optimization problem (20) as the MILP optimization problem (22) shown at the bottom of this page. This formulation is solvable by standard MILP solvers. The number of constraints grows as O(K + T ), as does the size of the optimization variables.
It is important to note that T is implicitly linked to p − λ max . Indeed, as we previously stated, we are looking for a solution that satisfies (20b). As we have assumed that λ n < p for all n ∈ {1, . . . , N}, then p − λ max ∈ [0, 1]; thus, there exists some z ∈ {0, . . . , T − 1} such that In the particular case where p is so small that 0 ≤ p − λ max < 1 T , then from (23) and the fact that |A x opt | ≤ T , we would have (24) In this case, the solver has to find a solution that accommodates every sample, and it will output the same solution regardless of the value of p. However, as one can see in Fig. 4, the rain attenuation CCDF, and subsequently the total attenuation, grows exponentially with log(1 − p). This exponential behavior has an impact on the size of the bandwidth used by the optimal carrier plan. Thus, having 0 ≤ p − λ max < 1 T is an issue, as the bandwidth required to accommodate all samples might be much bigger than what would be really needed for this value of p. In order to not oversize our carrier plan for a given p value, we then need to make sure that T ≥ 1 p−λ max .

B. Worst-Case Method
As a benchmark for our proposal, we consider a worstcase approach, which is simpler and does not consider the spatial correlation of fade. At worst, all terminals are always under fade simultaneously. Under this assumption, finding the optimal worst-case carrier plan x wc is simpler as the outage constraint is now only dictated by the link outage probability. The composition of x wc can be determined by finding, for each terminal n, the ModCod that will provide a link outage probability smaller than p. Traditional network sizing methods would achieve this by following the ITU recommendations and compute the attenuation CCDF for each terminal individually. Then, the CCDF is used to compute the corresponding SINR and, in turn, select the ModCod that will guarantee continuous operations at the fade levels corresponding to p. This method is commonly used in the satellite network dimensioning literature [7], [10]. In our method, we already computed time series of SINR. We can then define the empirical SINR CDF of each terminal n across time as for all n ∈ {1, . . . , N} and y ∈ R + in dB. The value of F SINR n (y) is the proportion of samples in which the SINR of terminal n falls below the value y. We can then compute C N+I wc n , the worst SINR met by terminal n with an outage probability p, with the generalized inverse distribution function (F SINR n ) −1 (p), that is We can then find the ModCod index K wc n ∈ {1, . . . , K} such that for all k ∈ {1, . . . , K}.

V. NUMERICAL RESULTS
In this section, we compare the MILP method, proposed in Section IV-A, with the worst-case method, described in Section IV-B. We first present the test scenario and parameters. We then compare the bandwidth produced by the carrier plans B milp and B wc , respectively, computed using the MILP and the worst-case methods.

A. Test Scenario
For our test scenario, we selected a Ka beam in Europe. We generated 489 rain maps of 200 by 200 km centered on (50.79 • N, 7.87 • E) with the MultiEXCELL model. In addition, N = 500 terrestrial terminals were generated in the same area based on the population density data of SEDAC [34] (see Fig. 5). For this article, we focused only on the rain component of the tropospheric attenuation as it is expected to be its dominant contribution. Fig. 4(a) and (b), respectively, shows the rain rate and the rain attenuation CCDFs obtained with such method for T = 30 000 Multi-EXCELL maps with solid lines representing averages and dotted lines representing minimums/maximums amongst terminals. On average, the synthetic curves are close to the ITU curves. The differences among the attenuation curves are due to the methods used to compute the attenuation from the rain rate. For the MultiEXCELL maps, we used the ATM PROP method [35], which differs slightly from the ITU Recommendation P.618 [36].
Other link parameters are given in Table II. As for the computation of the terminals' return link total SINR in (4), as we use a single beam, the term C IM n (t ) is expected to be relatively high. In addition, considering a star topology network, the downlink SNR C N dl n (t ) from the satellite to the gateway is also expected to be relatively high, assuming no fade on this link. On the other hand, C N ul n (t ) is expected to be relatively small due to limitations on the terminals' antenna size and amplifier power. Under those assumptions, the terminals' return link total SINR is, therefore, driven by the uplink SNR, such that We selected a maximum information rate of 1 Mb/s and a 1:50 contention ratio, as suggested in [37], giving us R CI = 200 kb/s. We used a ModCod set based on the DVB-RCS2 standard [38]. The ModCod SNR thresholds and spectral efficiencies, as given in Table III (K = 10), are extracted from [39, Table III].
Our results will show the evolution of both methods' performances across different values of p ∈ [p min = 1 · 10 −3 , p max = 5 · 10 −2 ]. As the sizing of EIRP 0 depends on the uplink rain margin, and thus on p, two different approaches were followed. 1) One where EIRP 0 is designed to ensure that the terminal with the most stringent link-budget in clear sky is able to afford the less demanding ModCod M 1 with an uplink rain margin equal to max n A n (p min ). Under this condition, the value of EIRP 0 , and consequently, λ max , does not depend on the actual value of p. Thus, the closer p is to p max , the bigger p − λ max becomes, and the terminals' transmission capabilities might be oversized. While this approach minimizes bandwidth costs, it maximizes the terminals' equipment cost (amplifier and antenna), and the total network cost might be suboptimal. We will refer to this approach as bandwidth over equipment (BOE). 2) A second one where EIRP 0 is recomputed at each given p to ensure that the terminal with the most stringent link-budget in clear sky is able to afford M 1 with an uplink rain margin equal to max n A n (p).
Under this condition, the values of EIRP 0 and λ max depend on the actual value of p, and the difference p − λ max stays (almost) the same across all values of p. This approach minimizes equipment costs but maximizes bandwidth costs. We will refer to this approach as equipment over bandwidth (EOB).
In a real network design, we expect that the bandwidth performances of our solution will fall in-between the results that are shown in the following sections.
In all scenarios, the value of λ max never exceeded 9.65 · 10 −4 . Thus, to satisfy the condition T ≥ 1 p−λ max , the 489 maps were circularly shifted multiple times in all directions across the area to obtain a number of time samples T = 30 000.

B. Bandwidth Performance
As depicted by the red curves in Fig. 6, the relative bandwidth gain of B milp w.r.t. B wc , i.e., (B wc − B milp )/B wc , increases when p decreases and becomes significant (greater than 10%) for p ≤ 1% in both scenarios. In both figures, the gain is negative for high p values (not commonly seen in VSAT networks). Then, it rapidly increases to positive values in the EOB scenario while the increase is slower in the BOE scenario. We then conclude that oversizing the uplink equipment will decrease the potential gain the MILP method offers over the worst-case method.
As many simulations had to be run, we imposed a time limit on the MILP solving. The consequence is that for some values of p, the solver was not able to find the optimal integer value. The black dotted curves in Fig. 6 represent the relative gap (B milp − B lb )/B milp between the bandwidth of the MILP solver's solution and the relaxed problem lower bound B lb . We observe that the gap curve is overall higher for the BOE scenario than for the EOB scenario. Thus, finding an actual minimum is likely to take more time when the uplink is oversized. However, a positive gap value does not necessarily imply that a better point exists but rather that not all options have been explored. Simulations for specific p values with longer time limits tend to show a reduction in gap due to finding a higher B lb rather than finding a more bandwidth efficient solution.
The microvariations in gain observed in Fig. 6 are mainly due to the variations of B wc , represented by a black  solid line in Fig. 7(a) and (b). While B wc follows a highly nonlinear trend, B milp shows an approximately linear behavior. We can then predict larger gains for outage probabilities lower than the ones showed. However, such probabilities are currently rarely used for the design of VSAT return links, and generating a solution would require a much larger number of samples.

C. Worst Terminal Outage
While the MILP method showed overall promising bandwidth performances, it is also important to verify the satisfaction of the outage probability constraint. To this end, we created a test algorithm that determines, for each samplê x(t ), which terminal will be dropped. For a given sample, the decision to drop terminal n is made by comparing C N+I n (t ) to C N+I wc n and the terminal's clear sky SINR. Terminals are dropped until the inequality in (13) is satisfied. Once x opt has been compared against all samples, the algorithm returns the highest outage proportion (link outage + drop due to congestion) amongst terminals. We will name this value the worst terminal outage p wto , and we expect p wto ≤ p. Fig. 8 shows the p wto obtained when testing the carrier plans produced by both methods. We can see that both the MILP and the worst-case methods are able to satisfy p wto ≤ p. Moreover, our MILP method provides lower p wto than the worst-case method. It is then apparent that increasing the bandwidth efficiency of the solution does not necessarily come at the expense of higher outage probability. The gap between the two methods is quite significant in the BOE scenario while moderate in the EOB scenario. We conclude that oversizing the uplink equipment has an impact on the oversizing of the MILP carrier plan.

A. MILP Solver Implementation
The gain of the MILP method depends on the solver's implementation and parameters. In our case, we used the integrated MATLAB integer linear programming solver intlinprog [40]. We limited the solver's execution in time to 6 h per value of p with simulations running in parallel on single core threads. We also instructed the solver to stop if it found an integer solution within 1% of the relaxed solution as gains lower or equal to 1% will not be significant in the context of network cost optimization. Execution times of several hours are affordable in our problem as we are doing long-term network planning. We also tested other solvers in specific scenarios, most notably SCIP [41] and CVX [42]. No notable differences in bandwidth gain or solver gap were observed with those parameters. However, we observed quicker reduction in the solver's gap when running a specific simulation on multiple cores with CVX.

B. Variability of Rain
In Fig. 6, we observed a crossing point at which the MILP bandwidth increased from negative values to positive ones. This crossing point depends on the average annual probability of rain p 0 across terminals. As p 0 depends on the terminals' geographical locations, we expect the crossing point to change for networks serving tropical or deserted regions, where p 0 is, respectively, higher and lower than in our European scenario.
The number of samples T also plays a role in the variability of rain amongst terminals. The lower the exceedance probability is, the lower the number of samples representing this probability is. As shown in Fig. 4, this translates into a high variability among terminals' CCDFs at low exceedance probabilities. Higher T values would reduce this effect and bring the dotted curves in Fig. 4 closer to the average curves.
The variability of rain amongst terminals is also impacted by the geographical separation between terminals in the network. In our scenario, we used 200 by 200 km maps as MultiEXCELL does not allow for the generation of much larger fields. The consequence is that the variability of rain is rather low and mainly tied to T . Expanding the network size would bring more variability and could change the performances of both methods. In future works, we plan to use different time series generators that span over larger areas, such as SISTAR.

C. Terminal Drop Policy
The curves shown in Fig. 8 depend on the terminal drop policy, of which an example was presented in Section V-C. One important property of this policy is that it stops dropping terminals when the accommodation property is satisfied. Thanks to the conservative nature of the accommodation property, the MILP method will guarantee p wto ≤ p for any dropping policy that satisfies the accommodation property. On the other hand, the worst-case method does not satisfy it; indeed, it can be straightforwardly verified by testing the resulting carrier plan on a drop policy which drops the terminals based on a predetermined order regardless of their SINR.

D. Heterogeneity for QoS and Number of Terminals
In a real VSAT network, it is common to have different QoS parameters (R CI and p) among terminals. In this context, a way to use our technique would be to group terminals under QoS groups, run the optimization for each group separately, and then merge the resulting carrier plans. This approach, however, has no optimality guarantees. Having an optimization problem able to deal with heterogeneous QoS requirements remains an open problem to be addressed in future works.
Another concern is the disparity in the number of terminals between the design stage and deployment phase, as well as its increase during the network lifetime. If the number of terminals was overestimated during the design phase, then the carrier plan is still valid, although not optimal. On the other hand, if the number of terminals was underestimated, then there will be more congestion than anticipated, and the QoS might not be achieved. As we already optimize the carrier plan for space-correlated fade, we expect our solution to be quite sensitive to an increase in the number of terminals. Therefore, real applications of our solution will likely have to anticipate a future rise in the number of terminals, and take margins accordingly.

VII. CONCLUSION
In this article, we presented the carrier plan sizing problem for adaptive VSAT return links and proposed a more bandwidth-efficient solution compared to current sizing methods. We defined the return link carrier plan design problem as applicable to broadband networks with adaptive return links and proposed an MILP optimization formulation that explicitly accounts for the spatial correlation of the tropospheric fade using time series generators. Finally, we introduced a sizing method commonly used in the literature and compared it to our MILP method in a realistic test scenario.
The numerical results illustrated how the proposed MILP method brings a bandwidth improvement between 10 and 50% for outage probabilities less than 1% while also being able to provide lower outage time than the benchmark method.
From those results, we concluded that the MILP method can improve the sizing of adaptive VSAT inbounds by reducing their cost and providing stronger QoS guarantees. Future works will focus on introducing hardware constraints, such as the terminals' power and antenna gain, and expand future test scenarios in scale and complexity (e.g., by introducing interference).