ConceptioArchivearXiv CS
arXiv CSopen access

Free-Placement Optimization of Ground Station Locations for Low-Earth Orbit Satellites

Unknown · 2026 · arxiv_cs
arXiv CS · Papers · License: Open Access · 2026
Open Source ↗Direct PDF ↓
distributedsystemsprotocols
networking, internet, protocols, distributed systems

Free-Placement Optimization of Ground Station Locations for Low-Earth Orbit Satellites

arXiv:2606.12667v1 [cs.NI] 10 Jun 2026

Grace Ra Kim∗ , Duncan Eddy† , Vedant Srinivas‡ , and Mykel J. Kochenderfer § Stanford University, Stanford, CA, 94305 Rapidly expanding low Earth orbit satellite constellations are placing increasing demands on terrestrial ground networks, motivating the development of more efficient ground station network designs. Current approaches select sites from predefined locations, limiting optimization to existing infrastructure and constraining performance. In contrast, free-placement optimization operates over a continuous spatial domain on Earth, broadening the search space and allowing higher-throughput configurations at the cost of potentially requiring new infrastructure deployment. In this work, we introduce SCORE (Sequential Cyclic Optimization via Refinement & Evaluation), a two-stage free-placement method for ground station design. SCORE combines sequential coordinate selection with cyclic refinement to manage high-dimensionality, nonconvexity, and local minima that challenge global optimizers. We benchmark SCORE against one-shot methods such as differential evolution (DE) and integer programming approaches using locations from Kongsberg Satellite Services and the World Teleport Association. Tests across two commercial Earth observation constellations (Capella Space and ICEYE) and one synthetic Walker-Star constellation show that SCORE requires up to 5× fewer function evaluations to converge relative to DE while improving downlink throughput by up to 13%. Compared to fixed-site methods, unconstrained SCORE achieves up to 15% greater total downlink, establishing a strong empirical performance benchmark for flexible placement; infrastructure-constrained SCORE retains over 92% of this gain while restricting placement to within proximity of existing fiber and power infrastructure. We also explore trade-offs between expanding existing stations and deploying new sites, informing future ground network design for operational constellations.

Nomenclature 𝑛 𝜆 𝜙 L 𝐿 S 𝑆 𝑓 𝑓data 𝑇opt start 𝑡 opt end 𝑡 opt 𝑇sim start 𝑡 sim end 𝑡 sim G

= = = = = = = = = = = = = = = =

number of selected stations longitude [degrees] latitude [degrees] set of ground station locations single ground station, coordinate (𝜆, 𝜙) set of satellites in a constellation single spacecraft in constellation performance objective function data downlink maximization objective function time period of mission duration [s] start time of mission duration [s] end time mission duration [s] time period of simulation [s] start time of simulation [s] end time simulation [s] set of all coordinate points on Earth

∗ Corresponding Author, PhD Student, Department of Aeronautics and Astronautics, Stanford, CA, AIAA Student Member, [email protected] † Postdoctoral Researcher, Department of Aeronautics and Astronautics, Stanford, CA, [email protected] ‡ Undergraduate Student, Department of Computer Science, Stanford University, Stanford, CA, [email protected] § Associate Professor, Department of Aeronautics and Astronautics, Stanford, CA, AIAA Associate Fellow, [email protected]

T I P 𝑃 𝑝 CL, S C𝐿 C𝑆 C𝑃 𝐶 𝑐 𝐿 dr 𝑆dr 𝐶dr 𝐶duration 𝐶 start 𝐶 end 𝑄 𝑄T 𝑄 proximity 𝑄I 𝐿T 𝑑min 𝑑 𝐿𝑖 ,𝐿 𝑗 𝑟 infra 𝐿I 𝑡 min 𝐸 max 𝑒 x 𝐷 𝑑 𝑁𝑃 𝐹 𝐶𝑅 u 𝐺 max 𝑔 X 𝑁𝜌 𝜌 Ω Ω𝜌 𝑠per plane 𝑠 𝑆 𝜌,𝑠 𝑀𝜌,𝑠 𝑓gap 𝑔gaps 𝑔¯ gaps 𝑛gaps

= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =

set of all coordinate points on Earth’s terrestrial surface set of feasible infrastructure on Earth’s surface set of predefined locations from existing Ground Station as a Service sites or teleport locations single ground station in P binary decision variable for selecting location 𝑃 contact windows between network L and constellation S set of contacts associated with station 𝐿 set of contacts associated with spacecraft 𝑆 set of contacts associated with location 𝑃 single contact object binary decision variable for selecting contact 𝐶 fixed ground station data rate [Gbps] fixed satellite data rate [Gbps] fixed contact data rate [Gbps] duration of contact 𝐶 [s] start time of contact [s] end time contact [s] penalty function penalty for land placement constraint penalty for station proximity constraint penalty for infrastructure proximity constraint single ground station site in T minimum distance between ground stations [km] distance between locations 𝐿 𝑖 and 𝐿 𝑗 [km] maximum allowable distance from infrastructure [km] closest point in I to station 𝐿 minimum contact duration for contacts [s] maximum number of cyclic evaluations in SCORE algorithm evaluation count in SCORE algorithm one-dimensional decomposed coordinate list of L for Differential Evolution algorithm dimensionality of x for Differential Evolution algorithm random index between 1 . . . 𝐷 for Differential Evolution algorithm population size parameter for Differential Evolution algorithm mutation scale factor for Differential Evolution algorithm crossover rate parameter for Differential Evolution algorithm mutated vector of x for Differential Evolution algorithm maximum generations for Differential Evolution algorithm generation count in Differential Evolution algorithm population of candidate solutions for Differential Evolution algorithm number of planes in Walker-Star constellation single plane in Walker-Star constellation Right Ascension of Ascending Node [degrees] Right Ascension of Ascending Node for plane 𝜌 [degrees] number of satellites in each plane of Walker-Star constellation satellite in-plane index of Walker-Star constellation Walker-Star satellite with plane and in-plane index 𝜌, 𝑠 mean anomaly of satellite 𝑆 𝜌,𝑠 in Walker-Star constellation [degrees] mean gap time minimization objective function all gaps in contact windows per-satellite mean gap time in contact windows count of all gaps in contact windows

2

I. Introduction The proliferation of satellite constellations in recent years is transforming global connectivity, planetary observation, and scientific exploration. Over the past decade, satellite deployment rates have accelerated markedly [1–3], driven by decreasing launch costs through ride-share missions [4] and the declining cost of small satellite development. This rapid growth has increased the demand for reliable ground station infrastructure that can meet latency, cost, and capacity requirements [5]. Selecting ground station locations poses a challenging decision for mission operators, as it is a multi-constrained problem driven by orbital alignment, access opportunities, cost, and site availability. Satisfying coverage, bandwidth, atmospheric conditions, infrastructure, and revisit requirements demands careful planning between satellite trajectories and ground segment assets, motivating the need for systematic location selection strategies. Despite its operational importance, ground station network optimization remains relatively underexplored in the literature. In this work, we solve the core ground station planning problem of selecting up to a predefined number of ground station locations that optimize the performance of a given mission objective. We categorize the ground station planning problem into two distinct variants: free placement and fixed-site selection. In the free placement problem, candidate locations may vary continuously in longitude and latitude, while the fixed-site selection problem restricts choices to a discrete set of predefined locations. These two variants address either end of the spectrum of planning environments, ranging from fully flexible deployments for newly designed networks to configurations constrained by existing infrastructure. We note that free placement considers only geometric (longitude and latitude) optimization, excluding terrain suitability (slope, obstructions) or civil engineering constraints characteristic of "greenfield" deployments. We focus particularly on ground station placement optimization in the context of serving low Earth orbit (LEO) Earth observation (EO) satellite constellations, which introduces the challenge of satellites continuously moving with respect to the station locations. This is in contrast to geostationary satellites, which are generally in continuous view of their supporting ground station network. The number of EO satellites in low Earth orbit has surged in recent years. As of 2024, there are 589 EO satellites in orbit, operated by 53 different commercial entities [6]. These satellites support a broad range of applications, including precision agriculture [7, 8], disaster monitoring for wildfire detection [9], and water resource management [10]. Due to the high-resolution sensors onboard, EO satellites produce significant volumes of imagery data. Most EO satellites produce over one terabyte (TB) of data per day [11], with some, such as NASA’s NISAR mission, delivering over 80 TBs daily [12]. These satellites typically operate in LEO between 500–650 km altitude, leading to short visibility windows with respect to individual ground stations [7]. Communication windows for LEO satellites typically last 3 to 10 minutes per access [13]. Collected data can be stored on board the spacecraft temporarily, but mission operators must maintain a sufficient downlink cadence to avoid buffer overflows, data loss, and data delivery delays that can severely impact operations [14]. Sub-optimal ground station placement can drastically reduce the number of these opportunities, further limiting communications capabilities and impacting overall mission performance. Prior work has largely studied optimization of spacecraft operations and contact scheduling given existing ground networks, rather than addressing the placement of ground infrastructure itself. Mixed Integer Linear Programming (MILP) is a commonly used technique in this space, offering optimal solutions to complex scheduling constraints. For example, the Shaving Model [15] used MILP to dynamically adjust pass durations across a fixed network, reducing missed connections by 23% over traditional integer programming approaches [16]. The China remote sensing satellite ground station applied both MILP and particle swarm optimization to break down multi-station relay tasks into tractable single-station sub-problems [17]. For constellations with 50 or more satellites, decentralized approaches such as geometric neighborhood decomposition have been introduced to localize decision-making while preserving global coordination [18]. Conflict resolution has also received significant attention: Terran Orbital, for example, implements a tiered priority system to allocate antenna time based on mission criticality [19], while others explore secure scheduling via quantum key distribution frameworks [20]. Despite these advancements in scheduling and operational efficiency, nearly all prior work assumes fixed ground station locations. As a result, the question of whether ground infrastructure is optimally placed to meet mission requirements remains largely unexamined. Ground station and gateway placement optimization literature commonly varies across three key dimensions: number of objectives, optimization techniques, and the problem formulation. Single-objective works predominate in fixed-site downselection, where meta-heuristics and heuristics select optimal subsets from predefined candidate pools. For instance, Guo et al. [21] employ particle swarm optimization, a meta-heuristic, to balance traffic load and deployment cost across fixed gateway candidates, while Liu et al. [22] use throughput evaluation heuristics for ground station planning from existing sites. Multi-objective formulations expand to grid-down selection and joint optimization challenges, often blending meta-heuristics with mathematical programming. Baeza et al. [23] integrate rain attenuation, visibility, elevation, and traffic in a multi-criteria grid-based framework for NGSO gateways; their follow-on work

3

[24] adds capacity diversity through similar grid downselection. Abe et al. [25] pursue joint placement-routing with mixed-integer programming on weighted multi-objectives of flow and cost, while Chen et al. [26] apply MILP heuristics to multi-gateway coverage in ISL-enabled constellations from candidate sets. Pure free-placement problem formulations, optimizations performed continuously over longitude and latitude without predefined sites, are less prevalent due to computational complexity. Architecture studies from Del Portillo et al. [27, 28] leverage genetic algorithms for optical/EHF ground segment design, trading availability, cost, and capacity while minimizing cloud cover and latency. A larger body of prior work has focused on selecting the best subset of stations from a predefined pool of existing candidate sites, rather than exploring methods for free placement of entirely new stations. In the EO context, Eddy et al. [29] formulate an integer programming (IP) based approach to identify data- and cost-optimal configurations using commercial Ground Station-as-a-Service (GSaaS) providers. For optical ground station (OGS) networks, Fuchs et al. [30] and others select the most effective subset of existing OGSs to support reliable space-to-ground optical communications [30, 31]. While candidate expansion and downselection suffice when infrastructure aligns with mission needs, terrain suitability, regulatory restrictions, power access, backhaul, and atmospheric factors often exclude optimal locations from the commercial pool. While pure continuous search methods address this gap, their computational cost hinders convergence for large networks; this necessitates scalable free placement optimization methods for greenfield networks supporting mission-scale LEO constellations. Outside of the space domain, a substantial body of literature exists on free placement location optimization in operations research [32], telecommunications [33], and urban management [34]. The optimization of warehouse, hospital, fire station, and cell tower locations to maximize coverage is a classic facility location problem, often addressed using meta-heuristic algorithms such as differential evolution [35, 36], genetic algorithms [33, 37], particle swarm optimization [38], and ant colony optimization [39]. These approaches navigate large, complex solution spaces, enabling flexible placement while balancing competing objectives such as cost, coverage, and demand. In such problems, gradient-based methods are often unreliable, as gradients may not exist at all feasible points or may fail to guide the search toward a global optimum. Similar characteristics appear in free placement ground station optimization problem, where discontinuous, non-convex objectives complicate the search. Adapting meta-heuristic or trial-and-error methods from other domains is therefore a natural approach, as they can explore these complex spaces and achieve near-optimal solutions [37, 40, 41]. Despite their effectiveness, even near-optimal solutions are computationally expensive in large search spaces, highlighting the need for more efficient methods [40, 42]. This work conducts a comprehensive study on free placement ground station optimization for LEO Earth observation satellite constellations. We introduce SCORE (Sequential Cyclic Optimization via Refinement & Evaluation), a framework for selecting ground station locations through iterative, cyclic coordinate search. Unlike traditional global optimization methods, SCORE decomposes the high-dimensional placement problem into a sequence of focused selection steps, enabling efficient exploration of large decision spaces. We evaluate SCORE’s performance through two comparative studies. First, we benchmark the algorithm against differential evolution, a free-placement metaheuristic method, where we demonstrate that SCORE achieves comparable or superior station placements while reducing the number of function evaluations needed for convergence by 5×. Second, to assess the advantages of continuous coordinate free-placement optimization over fixed-site selection, we compare SCORE to a globally optimal IP formulation that selects sites specifically from KSAT and locations listed by the World Teleport Association. SCORE’s continuous free-placement optimization identifies placement configurations that outperform fixed-site methods by up to 15%, highlighting the potential gains from flexible placement. Under operationally feasible latitude constraints, SCORE retains over 92% of this gain while remaining within commercially viable deployment regions. We demonstrate these advantages in case studies involving the 5-satellite Capella constellation and the 34-satellite ICEYE constellation, and explore trade-offs between deploying additional ground stations versus expanding antenna capacity at current infrastructure sites.

II. Problem Formulation The ground station placement problem is selecting a set of 𝑛 longitude–latitude (𝜆, 𝜙) coordinate locations L for a given satellite constellation S that maximizes a mission-specific performance objective 𝑓 . Ideally, the ground station placement problem could be posed as a scheduling optimization by predicting all possible satellite-to-ground contact opportunities over the entire mission duration, then identifying an optimal set of contact opportunities for this time period. This formulation would guarantee optimal performance over the entire mission duration, but due to the stochastic orbital perturbations experienced by LEO spacecraft, it is not possible to perform long-term orbital propagation with sufficient accuracy to predict all contact opportunities over multiple years.

4

To address this challenge, we adopt a surrogate optimization inspired by Eddy et al. [29] that focuses on performing ground station placement optimizations over shorter time periods, typically between 7 and 10 days, rather than the full multi-year-long mission duration. The cyclic nature of satellite orbits means that such a window captures multiple orbital cycles, producing a distribution of contact opportunities that is assumed to approximate the long-term contact opportunities of the constellation. While seasonal variations do affect the exact contact opportunities over long durations, the overall number and order over the shorter duration are assumed to be a representative surrogate. We perform a large-scale computational study in Appendix VI.C to validate this assumption. The full mission duration is defined by start and end time 𝑡 end , with total duration 𝑇 end start the optimization start time 𝑡opt opt = 𝑡 opt − 𝑡 opt . For surrogate optimization, opt start end a shorter simulation window is used, defined by start and end times 𝑡 sim and 𝑡 sim , respectively, resulting in a fixed end − 𝑡 start . one-week duration 𝑇sim = 𝑡sim sim We start our problem formulation by defining the set of all coordinate points G on Earth’s ellipsoidal surface, where  G = (𝜆, 𝜙) ∈ R2 | −180◦ ≤ 𝜆 ≤ 180◦ , −90◦ < 𝜙 < 90◦ ∪ {(0, −90), (0, 90)} (1) is the set of all possible longitude-latitude coordinate pairs (𝜆 𝑖 , 𝜙𝑖 ). A ground station network is then defined as L = {(𝜆1 , 𝜙1 ), . . . , (𝜆 𝑛 , 𝜙 𝑛 )},

L ⊆ G,

|L| = 𝑛

(2)

where L is the set of the final 𝑛 selected site locations. Every station 𝐿 ∈ L has a fixed data rate 𝐿 dr , representing the maximum number of bits that can be received per second at that site. In this work, 𝐿 dr is assumed constant across all sites, though it can generally vary between stations. The network selection is optimized with respect to a set of satellites S, where each spacecraft 𝑆 ∈ S has associated design constants such as its orbital elements (two-line-elements, or TLEs, in this work) and fixed data rate 𝑆dr . We study two problem formulations for selecting L, free placement and fixed-site selection. In the free placement problem, site locations L can be any longitude and latitude coordinate in G. This formulation models cases where ground station network locations need to be selected from any continuously varying location. Additional constraints can reduce the search space of G, such as limiting placement to land-based locations T where T ⊆ G, or further restricting to infrastructure-accessible zones I ⊆ T defined by proximity to existing population centers and terrestrial infrastructure. In this work, free placement optimization operates over geometric coordinate locations and does not account for terrain suitability or civil engineering constraints in traditional “greenfield” site selection. For the fixed-site selection problem, candidate ground station sites for the final ground network L are drawn from a discrete predefined coordinate list of ground stations P. Each location ( 𝑝, 𝑃) ∈ P has an associated binary decision variable 𝑝 ∈ {0, 1} that is 1 if the location is selected and 0 otherwise. Location 𝑃 is a coordinate (𝜆, 𝜙). The final set of selected ground station sites is then defined as L = {𝑃 | ( 𝑝, 𝑃) ∈ P, 𝑝 = 1}

(3)

the subset of all locations with 𝑝 = 1. This fixed-site selection problem is more applicable in the GSaaS or teleport selection cases, where existing ground station facilities or foundational infrastructures already exist and the selection problem is restricted to choosing from these preexisting locations. The selected networks are evaluated based on the quality of their possible contact windows CL,S , which represent all opportunities for communication between the network L and spacecraft S during the simulation window 𝑇sim . Since our goal is to estimate the achievable communication performance under realistic operational constraints, we formulate a scheduling problem over these contacts to approximate the network capacity achievable during the mission. Similar to P, each contact (𝑐, 𝐶) ∈ CL, S contains a binary decision variable 𝑐 ∈ {0, 1} that indicates whether that contact has been selected for scheduling. Each contact 𝐶 has a variety of other constants, such as the contact start time and end time (𝐶 start , 𝐶 end ). The data rate for each contact 𝐶dr is set as the minimum of the satellite and ground station data rates 𝐶dr = min(𝐿 dr , 𝑆dr )

(4)

participating in the contact. The duration of the contact is taken as the difference between the contact end and start times 𝐶duration = 𝐶 end − 𝐶 start .

5

A. Objective Functions Mission operators may prioritize different objectives when selecting a ground station network to support a satellite constellation. In general, these objectives can be expressed as selecting a set of sites that yields the most favorable contact windows CL, S between the spacecraft and the ground network, subject to the mission goals and constraints. We first define a general optimization formulation as max 𝑓 (L, S)

(5)

L⊂ G

where the decision variable is a candidate set of ground stations L ⊂ G and 𝑓 (L, S) is an objective function that evaluates the quality of the chosen network given the constellation S. This work focuses on designing ground station networks for high-data-volume EO constellations, with the primary goal of maximizing total downlinked data over the mission, independent of delivery latency. We call this the data downlink-maximization objective. To evaluate this objective for a given candidate network L, we first formulate a scheduling problem that determines which contacts between each spacecraft and the ground stations in L are feasible, subject to network and visibility constraints. Visibility is enforced through a minimum elevation angle requirement, such that a contact is feasible only if the satellite elevation exceeds 10◦ . Network constraints are further detailed in Sections II.B and II.C. Only contacts selected in the scheduling (𝑐 = 1) are included in computing the final objective. The total data downlinked from a network L is defined as 𝑓data (L, S) =

𝑇opt ∑︁ 𝑇sim 𝐿 ∈ L

∑︁

(6)

𝐶dr 𝐶duration 𝑐

(𝑐,𝐶 ) ∈ C𝐿 𝑇

opt where C𝐿 ⊆ CL, S represents the set of contacts associated with station 𝐿. We weight the objective by 𝑇sim to appropriately approximate the total data volume downlinked over the mission duration, rather than only the simulated time window. We also outline other objectives such as mean-gap minimization; details are provided in Appendix VI.A.

B. Regularization: Reducing the Search Space For the free placement problem, additional constraints need to be imposed on the search space G to prevent the selection of infeasible locations. We consider three constraints that define infeasible ground station placements: stations must be located on land; stations must maintain a minimum separation distance to prevent overlapping communication cones; and stations must be close to existing power and terrestrial communications infrastructure. However, if we incorporate every infeasible region as a hard constraint, we limit the optimization methods that can be used to solve the problem. Many efficient optimization techniques, such as gradient-based methods or global heuristics designed for unconstrained problems, are not naturally equipped to handle highly irregular, non-convex feasible regions, making it impossible or computationally expensive to enforce hard constraints directly. To address this, we incorporate these constraints as penalty terms subtracted from the objective function in Section II.A, discouraging violations without restricting the optimization method. Each penalty takes the form of a quadratic loss function 𝑄(𝐿) = (∥𝐿 − 𝐿 ∗ ∥) 2

(7)

where 𝐿 ∗ represents the nearest feasible ground station location that satisfies the constraint on 𝐿. When the constraint is satisfied, ∥𝐿 − 𝐿 ∗ ∥ = 0 and no penalty is applied. When violated, the penalty grows quadratically with the distance to the nearest feasible point. This provides a flexible tool for modeling ground station placement constraints, as any hard constraint can be incorporated into the objective function in this form. We first apply this to the land placement constraint. This constraint ensures that all ground stations can only be placed in feasible terrestrial locations L ⊆ T . Although this constraint focuses on excluding bodies of water, the definition of T can be modified to reduce the search space to exclude any infeasible station placement location. We transform this hard constraint into a penalty function like in Equation (7) as ∑︁ 𝑄 T (L) = (∥𝐿 − 𝐿 T ∥) 2 (8) 𝐿∈ L

where we calculate the geodesic distance between station 𝐿 and 𝐿 T , which represents the closest land mass coordinate to station 𝐿 in T . If station 𝐿 is in T , the distance between 𝐿 − 𝐿 T is zero. This penalty can be subtracted from our

6

objective functions in Section II.A, to penalize whenever a given ground station is not within the feasible terrestrial locations. The penalty grows quadratically the further station 𝐿 is from T . With this formulation of the quadratic loss penalty, any geographic region becomes heavily penalized in the search space, effectively discouraging its selection in the free placement ground station optimization problem. Next, we consider the station proximity constraint, which ensures a minimum distance 𝑑min is maintained between all ground stations so that network location redundancy is minimized and closely spaced ground stations do not experience interference limitations. In practice, this restriction may be relaxed through techniques such as polarization diversity or frequency allocation at high-demand ground sites, enabling multiple antennas to operate in close proximity. Incorporating such capabilities would require extending the formulation to jointly model spatial placement and communication resource allocation, which is beyond the scope of this work but represents an interesting direction for future research. Given a set of selected candidate locations L, 𝑑 𝐿𝑖 ,𝐿 𝑗 is denoted as the distance between locations 𝐿 𝑖 and 𝐿 𝑗 . for all locations 𝐿 𝑖 , 𝐿 𝑗 ∈ L. The value of 𝑑min is a predefined constant of the minimum allowable distance between any two stations. The constraint is formulated as 𝑑 𝐿𝑖 ,𝐿 𝑗 > 𝑑min

∀ 𝐿 𝑖 , 𝐿 𝑗 ∈ L where 𝐿 𝑖 ≠ 𝐿 𝑗

(9)

which ensures that any pair of candidate ground station locations 𝐿 𝑖 , 𝐿 𝑗 ∈ L, maintains a minimum separation distance of at least 𝑑min . To regularize this constraint, we frame Equation (9) again in the form of the general penalty function from Equation (7) 𝑛−1 ∑︁ 𝑛 ∑︁ 𝑄 proximity (L) = max(0, 𝑑min − 𝑑 𝐿𝑖 ,𝐿 𝑗 ) 2 (10) 𝑖=0 𝑗=𝑖+1

where 𝑄 proximity (L) is calculated between every combination of 𝐿 𝑖 , 𝐿 𝑗 ∈ L. The penalties associated with each 𝐿 𝑖 , 𝐿 𝑗 pair in Equation (10) can again be subtracted from the objective function 𝑓 (L, S) to apply these constraints. Finally, we consider an infrastructure proximity constraint, which encourages ground station placement near existing population centers and terrestrial infrastructure. This constraint is motivated by the practical requirement that ground stations require access to reliable power and fiber backhaul for high-throughput data downlink operations. We define a set of feasible infrastructure zones I ⊆ G as the union of circular regions of radius 𝑟 infra centered on known population centers, such that Ø I= {(𝜆, 𝜙) ∈ G | 𝑑 ((𝜆, 𝜙), (𝜆 𝑘 , 𝜙 𝑘 )) ≤ 𝑟 infra } (11) 𝑘

where (𝜆 𝑘 , 𝜙 𝑘 ) are the coordinates of known population centers and 𝑟 infra is the maximum allowable distance from infrastructure, set to 50 km in this work. Similar to the land placement constraint, we transform this into a penalty function as ∑︁ 𝑄 I (L) = (∥𝐿 − 𝐿 I ∥) 2 (12) 𝐿∈ L

where 𝐿 I represents the closest point in I to station 𝐿. The full regularized optimization problem from Equation (5) can now be written as max 𝑓 (L, S) − 𝑄 T (L) − 𝑄 proximity (L) − 𝑄 I (L) (13) L

with the regularization terms of the land placement, proximity between stations, and proximity to infrastructure penalty. For the fixed-site selection case, these penalties are set to zero as all locations in P are already pre-filtered to guarantee that they do not violate these conditions. With the above formulation, the free-placement problem is fully specified and can be solved using existing optimization techniques such as gradient-based search or evolutionary algorithms. However, for the fixed-site selection variation, additional work is needed to express the formulation as an integer program that explicitly models discrete decisions. C. Constraints and Fixed-Site Integer Programming Formulation After selecting objective and regularization terms, constraints are introduced to model system, design, and operational limitations. These can be grouped into two categories, contact exclusion constraints, which prevent infeasible overlaps in scheduling, and selection constraints, which enforce consistency between contacts, sites, and network size. These are necessary to properly compute the maximum feasible number of contacts that can be taken and data downlinked given a set of predicted contacts. In the free placement case, contact exclusion constraints are applied during evaluation after

7

candidate networks are selected, ensuring fairness when comparing to fixed-site results. In fixed-site selection, both contact exclusion and selection constraints are encoded directly into an integer program. This formulation provides a unified optimization framework: free placement problem formulations are optimized through penalty-regularized search in Equation (13), while fixed-site selection leverages integer programming formulation solvers with exact optimality guarantees. Together, they enable consistent evaluation of placement strategies under different operational contexts. 1. Contact Exclusion Constraints Contact exclusion constraints model the limitations of the spacecraft or ground station site’s ability to communicate with different assets. These constraints reflect the single-antenna limitation present either on satellites or at ground stations, ensuring that at any given time, each satellite communicates with at most one ground station and vice versa. To build the satellite contact exclusion constraint, all contacts for satellite 𝑆, denoted as C𝑆 , must be examined to identify any simultaneous contacts from multiple sites within candidate network L. For computational efficiency, the contacts can be sorted in ascending time-order to limit checks to temporally nearby contacts, though this is not required for the constraint itself. This constraint can be expressed as 𝑐𝑖 + 𝑐 𝑗 ≤ 1

∀ 𝑆 ∈ S, 𝑖, 𝑗 ∈ {1, . . . , |C𝑆 |}, 𝑗 > 𝑖 s.t.

𝐶𝑖start ≤ 𝐶 end 𝑗 end 𝐶 start ≤ 𝐶 𝑗 𝑖

(14)

where only one contact in every pair of overlapping contacts for a satellite is selected, leading to each satellite communicating with only one ground station at a time. The inverse of the satellite contact exclusion constraint can also be formulated to reflect ground station infrastructure limitations, where each station is restricted to serving only one satellite simultaneously. For the station contact exclusion constraint, the constraint is outlined as 𝑐𝑖 + 𝑐 𝑗 ≤ 1

∀ 𝐿 ∈ L, 𝑖, 𝑗 ∈ {1, . . . , |C𝐿 |}, 𝑗 > 𝑖 s.t.

𝐶𝑖start ≤ 𝐶 end 𝑗 start end 𝐶 𝑗 ≤ 𝐶𝑖

(15)

where each location is set to ensure communication can only occur with one satellite at any time. All contacts in C𝐿 are organized in time-order sequence. 2. Fixed-Site Integer Programming Selection Constraints The fixed-site ground station selection problem is particularly well suited for formulation as an integer program, as both the predefined ground station selection list P and set of contacts C(L, S) have associated binary decision variables indicating whether the contact (𝑐, 𝐶) or station site ( 𝑝, 𝑃) were selected. By coupling site-selection variables 𝑝 with contact-scheduling variables 𝑐, the IP formulation enables simultaneous optimization of station locations and communication schedules. This approach allows us to leverage established solver libraries such as Gurobi [43] or COIN-OR [44], which provide fast, efficient solutions to IP problems along with optimality certificates that indicate whether a solution is globally optimal, suboptimal, or infeasible. We follow existing integer programming selection formulations from prior works [29, 45], which solve the fixed-site ground station selection problem for GSaaS providers. To properly characterize the fixed-site selection ground station optimization problem using an IP, we consider the data downlink maximization objective defined in Section II.A and apply both contact exclusion constraints. We also introduce several additional IP selection constraints to ensure consistency between contact- and site-level decisions, enforcing network design requirements. The first additional constraint is to ensure that if a contact 𝑐 is selected, the decision variables of the location 𝑝 are also selected. No contacts should be scheduled at a location unless the location is activated. This is denoted as ∑︁ 𝑐 ≤ |C𝑃 | 𝑝, ∀ ( 𝑝, 𝑃) ∈ P (16) (𝑐,𝐶 ) ∈ C𝑃

where C𝑃 represents the set of contacts from site 𝑃. Second, we introduce a network size constraint to enforce the size

8

constraint 𝑛 on the ground station network ∑︁

𝑝 𝑖 = 𝑛, ∀ ( 𝑝, 𝑃) ∈ P

(17)

𝑖

to ensure we select ground station networks with the intended size. We also introduce a minimum contact duration constraint to ensure that only contacts with a duration greater than a specified 𝑡 min are considered for planning. The constraint is defined as (18) 𝑐 = 0 if 𝐶 end − 𝐶 start < 𝑡min , ∀ (𝑐, 𝐶) ∈ C which eliminates short-duration contacts that might not be operationally useful.

III. Methodology We outline methodologies for addressing the free placement ground station problem described in Section II. To navigate the large, high-dimensional search space, we employ two gradient-free optimization methods. First, we introduce Sequential Cyclic Optimization via Refinement & Evaluation (SCORE), which combines sequential cyclic search with single-step, gradient-free optimization methods; in this work, we demonstrate its application using Nelder-Mead [46] and Powell [47]. As SCORE is optimizer-agnostic, other gradient-free methods can be easily substituted for the single-step optimizer without changing the framework. Nelder-Mead was selected due to its simplicity, broad applicability in gradient-free settings, and well-understood convergence behavior in low-dimensional spaces. It is particularly suitable for local refinement, which aligns naturally with SCORE’s sequential cyclic updates. Powell was selected for its robustness to noisy or irregular objectives and its efficient line-search approach, providing a complementary local optimization strategy to Nelder-Mead within SCORE. Second, we review differential evolution (DE), a widely used evolutionary algorithm for location selection [35, 36]. DE is a population-based global optimizer that serves as a strong baseline for comparison. While it can achieve high-quality solutions, it requires substantially more objective evaluations than SCORE, especially as the number of stations grows. A. Sequential Cyclic Optimization via Refinement & Evaluation (SCORE) To address the computationally costly nature of one-step global optimization methods, we introduce SCORE, Sequential Cyclic Optimization via Refinement & Evaluation. SCORE decomposes the full optimization into a series of smaller problems and proceeds in two phases. (1) Sequential coordinate selection: beginning with an empty set, Algorithm 1 SCORE: Sequential Cyclic Optimization via Refinement & Evaluation 1: function SCORE(S, 𝑓 , 𝑛, 𝐸 max ) ⊲ 𝑛: desired # of ground stations, 𝐸 max : maximum # cyclic evaluations 2: Initialize L ← ∅, evaluation count 𝑒 ← 0 3: (1) Initial Sequential Coordinate Selection 4: while length of L < 𝑛 do 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19:

𝐿 ← Optimizer(L, S, 𝑓 ) ⊲ Optimizer can be any gradient free method L ← L ∪ {𝐿} ⊲ Add 𝐿 to set of ground station locations L (2) Sequential Cyclic Refinement while 𝑒 < 𝐸 max do Store copy of L as L1 for 𝐿 ∈ L do ⊲ New set L ′ without current station location 𝐿 L ′ ← L \ {𝐿} ′ ′ 𝐿 ← Optimizer(L , S, 𝑓 ) ⊲ Select new station 𝐿 ′ ′ ′ ′ ′ L ← L ∪ {𝐿 } ⊲ New set L with current station location 𝐿 if 𝑓 (L ′ , S) yields better objective value than 𝑓 (L, S) then L ← L′ if L1 = L then ⊲ No changes were made, current sequence of stations are optimal break 𝑒 = 𝑒+1 return L

9

SCORE sequentially adds ground stations to L, selecting at each step the longitude–latitude pair that maximizes the objective function given the current network configuration. This repeats until network L reaches the prescribed network size 𝑛. (2) Sequential cyclic refinement: with 𝑛 stations, the algorithm iteratively re-optimizes one station at a time while holding the others fixed, cycling through all stations to refine their positions and improve the objective, until a stopping criterion is met (e.g., negligible improvement or a budget on iterations). By reusing intermediate contact evaluations CL, S from modifying only one station 𝐿 ∈ L each time, SCORE reduces the number of contact evaluations of the entire network needed during the search process. This coordinate structure-aware strategy enables SCORE to retain solution quality while significantly lowering computational cost, making it more practical for large-scale or design-iteration-heavy scenarios. We define the full framework in Algorithm 1. In the first phase, sequential coordinate selection, the algorithm constructs the ground station set L incrementally. Starting from an empty set, it repeatedly adds a station location 𝐿 ∈ G that yields the largest improvement to the objective relative to the stations already chosen. This new location is added to L, and the process continues until the desired number of stations 𝑛 is reached. Each added location 𝐿 is locally optimal based on the existing set L, and the final combination of stations approaches near-optimal configurations at a fraction of the computation cost. SCORE’s second phase, sequential cyclic refinement, performs further fine-tuning of the combined set of stations. Here, we refine the station set through iterative updates. For each station 𝐿 ∈ L, the algorithm temporarily removes 𝐿 to form a reduced set L ′ = L \ {𝐿}, and re-optimizes that location using the same optimizer. The updated station 𝐿 ′ is then added back to form a new candidate set L ′ ∪ {𝐿 ′ }. The objective function values of the original set L and the modified set L ′ ∪ {𝐿 ′ } are compared, and if the new set improves performance, it replaces the current one. This process is repeated for each station 𝐿 ∈ L, and up to 𝐸 𝑚𝑎𝑥 total cycles. The optimization terminates early if a full pass through all stations results in no updates, indicating that the current set L has converged to a locally optimal solution. Any gradient-free optimizer can be used for the Optimizer noted for each coordinate selection in Algorithm 1. For this paper, we employ the Nelder-Mead simplex search algorithm [46] as well as Powell’s conjugate direction method [47]. Algorithm 2 Nelder-Mead (NM) 1: function NelderMead( 𝑓 , 𝛼, 𝛾, 𝜌, 𝜎, 𝑍 max ,S) 2: Initialize simplex {𝐿 1 , . . . , 𝐿 𝑘+1 } randomly; set 𝑧 ← 0 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19: 20: 21: 22: 23: 24: 25: 26:

while 𝑧 < 𝑍max and stopping condition not met do Sort simplex points {𝐿 1 , . . . , 𝐿 𝑘+1 } so that 𝑓 (𝐿 1 , S) ≤ · · · ≤ 𝑓 (𝐿 𝑘+1 , S) Compute centroid 𝐿 cen of best 𝑘 points (excluding worst) 𝐿 𝑟 ← 𝐿 cen + 𝛼(𝐿 cen − 𝐿 𝑘+1 ) if 𝑓 (𝐿 1 ) ≤ 𝑓 (𝐿 𝑟 ) < 𝑓 (𝐿 𝑘 ) then Replace 𝐿 𝑘+1 with 𝐿 𝑟 else if 𝑓 (𝐿 𝑟 ) < 𝑓 (𝐿 1 ) then 𝐿 𝑒 ← 𝐿 cen + 𝛾(𝐿 𝑟 − 𝐿 cen ) if 𝑓 (𝐿 𝑒 ) < 𝑓 (𝐿 𝑟 ) then Replace 𝐿 𝑘+1 with 𝐿 𝑒 else Replace 𝐿 𝑘+1 with 𝐿 𝑟 else if 𝑓 (𝐿 𝑟 ) ≥ 𝑓 (𝐿 𝑛 ) then if 𝑓 (𝐿 𝑟 ) < 𝑓 (𝐿 𝑘+1 ) then 𝐿 𝑐 ← 𝐿 cen + 𝜌(𝐿 𝑟 − 𝐿 cen ) else 𝐿 𝑐 ← 𝐿 cen − 𝜌(𝐿 cen − 𝐿 𝑘+1 ) if 𝑓 (𝐿 𝑐 ) < 𝑓 (𝐿 𝑘+1 ) then Replace 𝐿 𝑘+1 with 𝐿 𝑐 else for 𝑖 = 2 to 𝑘 + 1 of {𝐿 1 , . . . , 𝐿 𝑘+1 } do 𝐿 𝑖 ← 𝐿 1 + 𝜎(𝐿 𝑖 − 𝐿 1 ) 𝑧 ← 𝑧+1 return best 𝐿 ∈ {𝐿 1 , . . . , 𝐿 𝑘+1 } in simplex such that 𝑓 (𝐿, S) is minimized

10

⊲ (1) Reflection:

⊲ (2) Expansion:

⊲ (3) Outside contraction ⊲ (3) Inside contraction

⊲ (4) Shrink:

1. Nelder-Mead Similar to evolutionary algorithms, Nelder-Mead is a gradient-free direct search method with applications ranging from engineering [48], chemistry [49], control theory [50], and computer science [50]. The simplex, a polytope in 𝑘 dimensional space with 𝑘 + 1 points [51], is at the core of understanding how optimization is performed under Nelder-Mead. In the free-placement ground station problem, each vertex of the Nelder–Mead simplex corresponds to a candidate configuration of ground stations, denoted {𝐿 1 , . . . , 𝐿 𝑘+1 }, where each 𝐿 𝑖 in this simplex encodes the longitude and latitude of station 𝑖. Thus, a simplex {𝐿 1 , . . . , 𝐿 𝑘+1 } represents 𝑘 + 1 competing configurations of station placements in the search space. At each iteration, the simplex vertices are ordered by their objective values, and the worst point is reflected across the centroid 𝐿 cen of the remaining 𝑘 points to explore a potentially better location. The behavior of this transformation is governed by four parameters: the reflection coefficient 𝛼 > 0, expansion coefficient 𝛾 > 1, contraction coefficient 𝜌 ∈ (0, 1), and shrinkage coefficient 𝜎 ∈ (0, 1). If the reflection improves the objective sufficiently, the new point is accepted; otherwise, the algorithm attempts an expansion (if the reflection is very good) or a contraction. Two types of contraction exist: outside contraction, when the reflection is moderately good, and inside contraction, when the reflection is worse than the current worst. If neither contraction yields improvement, the algorithm shrinks the entire simplex toward the best point, scaled by 𝜎. This iterative process continues until convergence or a maximum number of iterations 𝑍max is reached. The best vertex is then returned as the approximate solution. We provide the full pseudocode of Nelder-Mead in Algorithm 2. We make two implementation adjustments to the Nelder-Mead algorithm to appropriately adjust the method for use in the ground station placement problem. First, rather than optimizing in a two-dimensional longitude–latitude space (𝑘 = 2), we project candidate locations onto the three-dimensional unit sphere and perform the optimization in that domain (𝑘 = 3). This requires a simplex with 𝑘 + 1 = 4 vertices (a tetrahedron) instead of a triangle. Operating directly in three dimensions avoids degeneracies at the poles that arise in the longitude–latitude system, where all longitudes collapse to a single point at 𝜙 = ±90°. Second, we select randomized starting simplexes on the unit coordinate sphere with the largest possible simplex area. Specifically, we generate 100 random feasible points on the unit sphere and evaluate all (𝑘 + 1)-point combinations to form candidate simplexes. The simplex with the largest volume is selected, ensuring a well-scaled and well-distributed starting configuration. Without such a spread, Nelder-Mead can become trapped in poor local minima due to insufficient coverage of the solution space. 2. Powell’s Method Powell’s conjugate direction method [47] is a derivative-free direct search optimizer that minimizes a function of 𝑚 variables through successive one-dimensional line minimizations. Similar to Nelder-Mead, Powell is a gradient-free method that only requires function evaluations. The method maintains a set of 𝑚 search directions r1 , . . . , r𝑚 , initialized to the standard basis vectors e1 , . . . , e𝑚 . At each outer iteration, the algorithm sequentially minimizes 𝑓 along each direction in turn, updating the current point w after each line minimization. The net displacement v = w − p accumulated across all 𝑚 line searches is then used as a final line minimization step, exploiting the aggregate direction of improvement across the full cycle. The key innovation of Powell’s method is its direction update rule, which progressively builds a set of mutually conjugate directions. After each full cycle, the oldest direction r1 is discarded, the remaining directions are v shifted, and the normalized displacement |v| is appended as the new direction r𝑚 . For a strictly quadratic objective, this procedure guarantees that after 𝑚 complete cycles, all search directions are mutually conjugate with respect to the Hessian, and the exact minimum is found [47]. For non-quadratic objectives such as 𝑓 in this work, the method provides fast practical convergence without requiring derivatives. The algorithm terminates when the displacement norm |v| falls below a tolerance 𝜀 pow or a maximum number of iterations 𝑍max is reached, returning the best point w found. We provide the full pseudocode in Algorithm 3. B. Differential Evolution Differential evolution [52] is a robust, non-gradient-based meta-heuristic capable of efficiently exploring complex, high-dimensional solution spaces, making it suitable for a wide range of global optimization problems. The method is similar to evolutionary approaches such as genetic algorithms, but DE uses a population of real-valued vectors rather than encoded genotypes. With a real-valued vector population, the mutation and recombination steps can occur directly in continuous solution spaces without requiring encoding or decoding. New vectors in the population are generated from combinations of existing candidates, and the DE algorithm always retains the best-scoring solution. This black-box approach does not require gradient information to guide optimization. To apply DE to the free placement ground station problem, the ground station candidate list L is decomposed

11

Algorithm 3 Powell’s Conjugate Direction Method 1: function Powell( 𝑓 , 𝑍 max , 𝜀 pow , S) 2: Initialize starting point w; set search directions r𝑖 ← e𝑖 for 𝑖 = 1, . . . , 𝑚; set 𝑧 ← 0 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15:

while 𝑧 < 𝑍max and stopping condition not met do p←w ⊲ Save starting point for this cycle for 𝑖 = 1 to 𝑚 do 𝑡 ∗ ← arg min𝑡 𝑓 (w + 𝑡 r𝑖 , S) ⊲ Line minimization along r𝑖 w ← w + 𝑡 ∗ r𝑖 v←w−p ⊲ Net displacement (new conjugate direction candidate) if ∥v∥ < 𝜀pow then break ⊲ Convergence ∗ 𝑡 ← arg min𝑡 𝑓 (w + 𝑡 v, S) ⊲ Line minimization along v w ← w + 𝑡∗ v Update directions: r𝑖 ← r𝑖+1 for 𝑖 = 1, . . . , 𝑚 − 1; r𝑚 ← v/∥v∥ ⊲ Drop r1 , append v 𝑧 ← 𝑧+1 return w such that 𝑓 (w, S) is minimized

into a one-dimensional vector of ordered longitude and latitude coordinates, denoted as x. The algorithm maintains a 𝑁𝑃 population of such candidate solutions X = {x𝑖 }𝑖=1 , of population size 𝑁 𝑃, balancing exploration with computational cost. New candidate vectors are generated using the mutation scale factor 𝐹 ∈ [0, 2], which sets the magnitude of modifications to each vector, governing the trade-off between global exploration and local refinement at each iteration. Components of these mutant vectors may replace the current candidates according to the crossover rate 𝐶 𝑅 ∈ [0, 1], where 𝐶 𝑅 = 1 means every component of x is replaced, and 𝐶 𝑅 = 0 means no replacements occur. The process continues for a maximum number of generations 𝐺 max , providing a practical stopping criterion. DE requires access to the objective function 𝑓 , the dimensionality 𝐷 of x (equal to 2𝑛, as each of the 𝑛 ground stations is represented by a longitude and latitude coordinate), and the bounds for each dimension in x. The full pseudocode is provided in Algorithm 4. Algorithm 4 Differential Evolution (DE) 1: function DifferentialEvolution( 𝑓 , 𝑁 𝑃, 𝐹, 𝐶 𝑅, 𝐺 max , 𝐷, S) 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19: 20: 21:

𝑁𝑃 Initialize population X = {x𝑖 }𝑖=1 of dimensionality 𝐷 randomly within bounds Initialize generation count 𝑔 ← 0 while 𝑔 < 𝐺 max and solution not converged do for each base vector x𝑖 in population do (1) Mutation: Randomly select 3 distinct vectors such that x1 ≠ x2 ≠ x3 ≠ x𝑖 v𝑖 = x1 + 𝐹 · (x2 − x3 ) (2) Crossover: Pick a random vector index 𝑑 ∈ 1, . . . , 𝐷 for 𝑗 = 1, . . . , 𝑑 do Generate rand 𝑗 ∼ 𝑈 (0, 1) if rand 𝑗 < 𝐶 𝑅 or 𝑗 = 𝑑 then u𝑖, 𝑗 = v𝑖, 𝑗 else u𝑖, 𝑗 = x𝑖, 𝑗

(3) Selection: if 𝑓 (u𝑖 , S) ≤ 𝑓 (x𝑖 , S) then x𝑖 ← u𝑖 Update generation count: 𝑔 ← 𝑔 + 1 return best x𝑖 in population such that 𝑓 (x𝑖 , S) is minimized

12

The general DE outline consists of performing a mutation, crossover, and selection process for each vector x𝑖 in X. The mutation is computed using three distinct vectors in the maintained population, separate from the vector x𝑖 being modified. The scaling factor of this mutation is dictated by 𝐹. Each dimension of vector x𝑖 then undergoes the crossover step, where each index either keeps their original value, or crosses over with the mutated vector value based on the crossover rate 𝐶 𝑅. A random vector index 𝑑 ∈ 1, . . . , 𝐷 is selected beforehand to ensure that in the rare case that all indexes randomly are not selected to cross over, the index 𝑑 in vector x𝑖 is ensured to have a mutation. The mutated vector is saved in u𝑖 . Finally, the algorithm selects between the original vector x𝑖 and the modified vector u𝑖 , selecting the vector with the lower objective value. If the new vector has a more optimal objective function value, this vector replaces the existing vector x𝑖 in the population. After the solution reaches convergence, or the maximum number of generations 𝐺 max is completed, the candidate vector x𝑖 with the lowest objective value is returned. Although DE’s non-gradient approach allows for the search of the discontinuous search space of ground station placement optimization, the method’s reliance on a large number of function evaluations makes the method computationally expensive. For problems where each computation evaluation is lightweight, this overhead may be acceptable. However, in our case, where every candidate configuration x𝑖 , or ground station 𝐿 requires recomputing contacts C𝐿,S for evaluation, the computational cost quickly becomes a bottleneck. In addition, the number of candidate solutions needed to be evaluated grows with the size of the ground station network L needed to be optimized, which limits DE’s scalability for larger network design problems. SCORE addresses these computational challenges, as a more efficient, gradient-free alternative tailored to the structure of the ground station placement problem.

IV. Experimental Setup Before discussing the performance of SCORE against other ground station optimization methods, we outline our experimental simulation parameters for reproducibility. For all experiments, we used the simulation parameters outlined in Table 1 and propagated each spacecraft’s dynamics using the SGP4 propagator [53]. Following the surrogate optimization framework in Section II, we restricted the simulation to the first 7 days (𝑇𝑠𝑖𝑚 ), which was scaled by 𝑇𝑜 𝑝𝑡 ◦ 𝑇𝑠𝑖𝑚 to approximate the full mission duration 𝑇𝑜 𝑝𝑡 . We applied a 10 elevation mask at all ground station locations when calculating contact opportunities. In the current formulation, we assume a fixed data rate for all ground stations and satellites throughout each contact window, independent of elevation angle, and do not model frequency band assignments. All simulations were run on a dual-socket workstation with a total of 28 physical cores (56 threads), using 2.6 GHz Intel Xeon E5-2690 v4 processors and 128 GB of RAM. Table 1

Simulation Parameters

Simulation Parameters 𝑠𝑡 𝑎𝑟𝑡 𝑡 𝑠𝑖𝑚 𝑒𝑛𝑑 𝑡 𝑠𝑖𝑚 𝑡 𝑜𝑠𝑡𝑝𝑡𝑎𝑟𝑡 𝑡 𝑜𝑒𝑛𝑑 𝑝𝑡 𝐶𝑑𝑟

Values 2025-04-01 17:23:40.69 UTC 2025-04-08 17:23:40.69 UTC 2025-04-01 17:23:40.69 UTC 2026-04-01 17:23:40.69 UTC 1.2 Gbps

A. Satellite Constellations In our evaluations, we ran experiments on both synthetic constellations and real satellites, with orbits defined by Two-Line Element (TLE) sets. We used synthetic constellations specifically in parameter sweeps that explored the impact of varying constellation sizes. By generating synthetic constellations, we could precisely adjust simulation environments to measure how increasing constellation size affected the optimizer’s speed to reach convergence. For experiments that did not involve varying constellation sizes, we used existing EO satellite constellations, with TLEs providing orbital data for active satellites. In practice, satellite operators select ground stations from a limited set of locations offered by GSaaS providers to establish a supporting ground network. In this context, we compared the performance of our unconstrained free-placement approach with fixed-site selection, which represents current industry practice. By directly evaluating fixed-site selection, which employs IP-based solvers constrained to existing ground station locations, against the free-placement approach, we can quantify the potential gains achievable when the 13

limitations of pre-defined site locations are removed. 1. Synthetic Walker-Star Constellations In the unconstrained free-placement problem evaluation, we performed optimizations over a Walker-Star constellation with the parameters listed in Table 2. This configuration was used to perform parameter sweeps over the total number of satellites in the simulation, ranging from 1 to 4. For our Walker-Star formation, we restricted each orbital plane to contain only one satellite, so the total number of orbital planes 𝑁𝜌 equals the number of satellites |S|. Table 2

Walker-Star Constellation Parameters

Walker-Star Constellation Parameters

Values

Altitude Eccentricity Inclination Number of Planes Varied in Parameter Sweep (𝑁𝜌 ) Satellites per Plane (𝑠per plane )

781 km 0.001 86.4° 1–4 1

2. Existing Earth Observation Constellations For the existing EO satellite constellations considered in this paper, we focused on two major EO operators: Capella Space and ICEYE. The associated orbital characteristics are listed in Table 3. We chose these two constellations specifically for the variance in their orbital parameters. Capella Space represents a smaller, mixed constellation of both mid-inclination (45◦ , 53◦ ) and sun-synchronous (97◦ ) orbits, while ICEYE represents a more traditional sun-synchronous-only constellation of a larger size. Examining these two constellations enables us to identify the distinct challenges that arise when designing ground station networks for mixed-orbit versus sun-synchronous trajectories. Table 3 EO Satellite Operators Capella Space ICEYE

Earth Observation Satellite Constellations Constellation Size

Orbital Altitudes

Inclination Angles

5 34

525-575 km 560-580 km

45◦ , 53◦ , 97◦ 97◦

B. Existing Ground Station Networks To compare SCORE against methods that optimize over fixed-site selections, we considered two datasets of ground station locations as potential fixed sites, depicted in Figures 1 and 2. The first dataset is provided by Kongsberg Satellite Services (KSAT), the largest currently existing GSaaS provider. Selecting sites from KSAT’s extensive network represents a particularly challenging GSaaS optimization problem. We also considered a geographically diverse subset of 100 locations from the World Teleport Association’s (WTA) global teleport map. The full WTA dataset of 915 stations exceeds the scale at which exact IP solutions can be obtained within a reasonable time, and we therefore selected 100 stations to maximize geographic diversity while maintaining computational tractability. This dataset enables evaluation of network planning in regions with existing satellite communication infrastructure, where operators may wish to expand or optimize ground station deployments even in the absence of a dedicated GSaaS provider. Scaling fixed-site selection to larger candidate sets remains an open challenge, though recent works have begun to address scalable ground station network optimization [45].

V. Results We organize the results into two main evaluations. First, we compare SCORE with differential evolution to assess computational efficiency, convergence, and the potential for multiple near-optimal solutions using a synthetic Walker-Star 14

90 60 30 0 −30

0

0

0

12

15

18

90

60

30

0

0 −3

0 −6

0 −9

20 −1

−1

−1

80

−90

50

−60

Fig. 1 KSAT Ground Station Network. Locations of operational ground stations from KSAT. Communications cones assume 525km altitude and 10° minimum elevation angle.

constellation. This analysis quantifies SCORE’s convergence speed improvements relative to DE, showing that SCORE achieves equivalent or superior objective values while reducing the number of function evaluations to convergence by up to 5× fewer function evaluations than DE for larger ground station selection configurations. Next, we evaluate SCORE against fixed-site selection methods based on integer programming (IP) to analyze performance tradeoffs when network locations are constrained to existing infrastructure. Using the KSAT and WTA teleport networks as baselines, we quantify how limited site availability restricts achievable downlink performance compared to SCORE’s unconstrained placement. Focusing on site placement for the CAPELLA and ICEYE constellations, we found that SCORE consistently outperformed IP-based fixed-site methods, achieving improvements in data throughput of over 9% relative to fixed-site solutions. Collectively, these evaluations highlight SCORE’s scalability, efficiency, and applicability across a variety of ground station network design scenarios. A. Free Placement Problem We begin by comparing methods for the free-placement problem, benchmarking SCORE against differential evolution using the synthetic Walker-Star constellation described in Section IV.A.1. The mission objective 𝑓 was to maximize total data downlink over the simulation period, 𝑇𝑠𝑖𝑚 as outlined in Equation (6). Final objective function 𝑇 𝑝𝑡 values were scaled by 𝑇𝑜𝑠𝑖𝑚 to approximate the data downlink over entire optimization period. We varied the number of satellites and ground stations independently over {1, 2, 3, 4}, generating 16 unique scenarios to evaluate algorithmic performance and computational efficiency. Our evaluations focused on the number of function evaluations to reach convergence and solution quality, with each scenario repeated over 10 random seeds to ensure robustness. Contact exclusion constraints were handled in two steps to prevent overlapping satellite contacts across ground stations. After optimization, solutions were post-processed with an integer programming solver to enforce strict non-overlapping assignments, counting any overlaps only once in the final objective. 1. Computation and Performance Gains First, we evaluated SCORE and DE’s tradeoffs between computation efficiency and mission objective performance. Figure 3 displays two critical metrics: (a) the number of function evaluations required for each algorithm to reach solution convergence and (b) the final objective values of total data downlinked in petabytes over 𝑇𝑜 𝑝𝑡 . We evaluate SCORE using two derivative-free local search methods, Nelder-Mead and Powell, demonstrating optimizer flexibility

15

90 60 30 0 −30

18 0

15 0

12 0

90

60

30

0

0 −3

0 −6

0 −9

20 −1

−1

−1

80

−90

50

−60

(a) Global Teleport List. Locations of a geographically diverse subset of 100 teleport stations from the World Teleport Association’s global teleport map, selected to maintain computational tractability of the integer programming formulation. Communications cones assume 525km altitude and 10° minimum elevation angle.

90 60 30 0 −30

0 18

0 15

0 12

90

60

30

0 −3

0

0

20 −1

−6

50 −1

0

80

−1

−90

−9

−60

(b) Full Global Teleport List from the World Teleport Association’s global teleport map. Communication cones are not plotted for better visibility of locations.

Fig. 2

Global Teleport List Maps.

16

SCORE (Powell)

Differential Evolution

0.98

1

1

1

0.46

0.47

0.45

0.48

0.93

0.87

0.78

0.74

2

2.1

2.1

2

0.95

0.88

0.89

0.87

2.8

3.2

2.4

2.9

3

3.1

3

3.1

1.4

1.4

1.4

1.3

10

11

10

5.7

3.9

4.1

4

4.1

1.9

1.8

1.8

1.8

16

17

21

17

2 3 4 Number of Satellites

1

2 3 4 Number of Satellites

1

10

1

1

2 3 4 Number of Satellites

# of Function Evaluations (103, log)

Number of Ground Stations 4 3 2 1

SCORE (Nelder-Mead)

(a) Total number of function evaluations to reach convergence for SCORE (Nelder-Mead and Powell) and DE free placement algorithms SCORE (Powell)

Differential Evolution

0.5

1.4

1.0

0.1

0.3

0.6

0.7

0.5

0.5

1.4

1.0

0.8

1.0

2.5

1.9

0.3

0.5

0.8

1.0

0.9

1.0

2.6

2.0

3.0 2.5 2.0

0.9

1.3

2.9

2.6

0.3

0.9

1.5

1.7

0.9

1.2

2.6

2.7

0.9

1.6

3.3

3.4

0.5

1.0

1.7

2.1

0.8

1.5

3.2

3.2

2 3 4 Number of Satellites

1

2 3 4 Number of Satellites

1

1.5 1.0

1

0.5

Data Downlinked in PB over Topt

Number of Ground Stations 4 3 2 1

SCORE (Nelder-Mead) 0.4

2 3 4 Number of Satellites

(b) Final data downlinked in PB for SCORE (Nelder-Mead and Powell) and DE free placement methods

Fig. 3

Computation versus performance gains for SCORE versus leading free placement optimization algorithms

while maintaining lower per-iteration cost than global meta-heuristic methods such as DE. In the heatmaps of Figure 3a, SCORE’s number of function evaluations increases linearly with the number of ground stations 𝑛 due to its sequential selection process. However, DE’s number of function evaluations grows much more rapidly, reaching over 17,000 function evaluations in the largest scenario compared to SCORE’s 4,100. For both algorithms, varying the constellation size |S| had minimal effect on number of function evaluations needed until convergence. Figure 3b shows the final mission objective values of total data downlinked in petabytes for both methods. Crucially, SCORE’s efficacy is sensitive to the choice of underlying optimizer; while the Powell method offers the lowest computational overhead (averaging roughly 50% fewer evaluations than Nelder-Mead), it results in significantly degraded objective values, often failing to surpass DE performance on even the simplest of configurations. Conversely, using a different optimizer like Nelder-Mead allows SCORE to maintain its predictable linear growth in complexity while consistently attaining near-optimal downlink performance. In larger networks, SCORE with Nelder-Mead surpasses DE by as much as 13% in total data downlinked (1 satellite, 4 ground stations), despite DE requiring up to five times more function evaluations (e.g., 21,000 vs 4,000 evaluations at 4 ground stations and 3 satellites). These results highlight SCORE’s superior efficiency, scalability, and effectiveness in exploring the solution space for ground station network optimization, provided a high-performing local search method is employed. 2. Convergence Behavior To better understand differences in convergence behavior between SCORE and DE, Figure 4 shows performance over time for the largest scenario from Figure 3: a four-satellite Walker-Star constellation requiring placement of four ground stations. We plot the total data downlinked against the number of function evaluations completed by each method. SCORE quickly reaches near-maximum downlink, achieving a mean performance of 3.4 PB in an average of 4,160 function evaluations (orange curve, vertical dashed line). DE’s increase is slower and more variable, reaching only 3.2 PB after an average of 17,432 function evaluations (blue dashed line). SCORE’s fast convergence can be attributed to its strategic sequential coordinate selection phase and Nelder-Mead’s efficient simplex search, which require few function evaluations to iteratively improve solutions. After fast selection of

17

2

17,432

4,160.0

Data Downlinked in Topt (PB)

3

SCORE Min/Max Range

1

DE Min/Max Range SCORE Final Evals DE Final Evals

0 0

2500

5000 7500 10000 12500 Number of Function Evaluations

15000

17500

Fig. 4 Data downlink comparison for SCORE and DE. SCORE (orange) converges towards a solution in around 4000 function evaluations; DE (blue) converges more slowly with higher variability.

an initial set of coordinates, the subsequent cyclic refinement facilitates incremental improvements until convergence, finding solutions with up to 5× fewer function evaluations than DE, while achieving comparable solution quality. The shaded regions surrounding each algorithm’s average convergence curve reflect variability across random seeds. SCORE’s variability remains limited, while DE’s stochastic, population-based search manifests in broader variability and less predictable convergence rates, due to its reliance on initial populations and random genetic operations. For further ablation studies, we direct the reader to Appendix VI.D, where we explore a wider variety of DE parameters, and Appendix VI.E, where we validate the cyclic refinement step and its robustness to selection ordering, both for the largest scenario of 4 ground stations and 4 satellites. 3. Multiple Near-Optimal Solutions In experiments varying random seeds, we found multiple distinct ground station layouts for the four-satellite Walker-Star constellation with four stations. Several example arrangements are displayed in Figure 5. Despite differences in spatial arrangement, all configurations achieved similar objective values, indicating multiple local optima. This suggests that overall layout patterns, rather than exact station positions, are the key factors influencing performance. The full list of coordinates and their data downlink values over 𝑇𝑜 𝑝𝑡 is provided in Table 4. Table 4 Aligned ground station coordinates from three optimization cases, reordered to highlight spatial similarity. Overall distribution patterns, not exact sites, primarily drive data downlink performance.

Scenario A B C

Coordinate 1

Coordinate 2

Coordinate 3

Lon.

Lat.

Lon.

Lat.

Lon.

Lat.

Lon.

Lat.

Data (PB)

15.65 −26.51 25.75

78.23 64.14 71.17

2.53 2.53 2.53

−72.01 −72.01 −72.01

−133.72 −148.49 −51.72

68.36 70.26 64.18

−57.85 168.38 −70.87

−51.68 −46.53 −52.94

3.11 3.10 3.12

18

Coordinate 4

90 60 30 0 −30

0

0

18

0

15

12

90

60

30

20 −9 0 −6 0 −3 0 0

50

−1

−1

80

−90

−1

−60

−60

−90

−90

−1

−1

−1

−1

−1

−1

(b)

12 0 15 0 18 0

−60

90

−30

60

−30

30

0

80 50

0

12 0 15 0 18 0

30

90

30

60

60

30

60

20 −9 0 −6 0 −3 0 0

90

80 50

90

20 −9 0 −6 0 −3 0 0

(a)

(c)

Fig. 5 Multiple high-performing ground station layouts yield similar objectives, showing placement flexibility. Coordinates and objective values are listed in Table 4.

19

B. Free Placement versus Fixed-Site Selection Methods We next evaluate SCORE’s performance against fixed-site selection optimization using integer programming solvers. To introduce greater realism, we compare SCORE to IP optimizations over existing ground station networks, specifically those provided by KSAT and WTA’s global teleport map. We evaluate three variants of SCORE against our fixed-site methods. First, unconstrained SCORE performs free placement with no infrastructure penalties, serving as a strong empirical performance benchmark on achievable downlink performance. Second, SCORE (Lat-Const) restricts placement to latitudes within [−73◦ , 78◦ ], corresponding to the range of existing KSAT infrastructure, representing a middle ground between fully unconstrained and infrastructure-constrained optimization. Third, SCORE (Infra-Const) incorporates the infrastructure proximity penalty from Equation (12), encouraging placement near existing population centers with access to fiber and power infrastructure. Comparing these three variants isolates the contribution of purely geometric optimization from the practical costs of site feasibility, quantifying the performance tradeoff operators face when moving from idealized to operationally realistic deployments. We further focus our evaluation on maximizing data downlink for established Earth observation satellite constellations, namely Capella Space and ICEYE, with their characteristics detailed in Section IV.A.2. When optimizing over a predefined set of candidate locations P, the IP solver returns a globally optimal solution, assuming a feasible solution exists. 1. Performance between SCORE, Fixed-Site Selection from Teleports and KSAT The downlink performance of SCORE and IP-based optimization over the mission window 𝑇𝑜 𝑝𝑡 is shown in Figure 6 for networks with 1 to 20 stations. For the Capella Space constellation (Figure 6a), unconstrained SCORE consistently outperforms IP-based methods, although IP achieves similar results for smaller networks. The strong performance of the fixed-site IP solutions can be attributed to Capella’s mid-latitude orbit inclinations, which align well with the predominantly mid-latitude distribution of ground stations and teleports shown in Figures 1 and 2. In this case, unconstrained SCORE yields only marginal gains for small networks, 0.4% over teleports and 4.6% over KSAT for three stations, but its advantage grows with network size, reaching improvements of 7% over teleports and 22% over KSAT at 20 stations. For the ICEYE constellation (Figure 6b), unconstrained SCORE significantly outperforms both IP baselines, leading by 8–15% across all network sizes. This reflects the limited high-latitude ground stations in KSAT and teleport networks, which hinders support for ICEYE’s high-inclination orbits. We note that these unconstrained results represent a strong empirical performance benchmark, as this version of SCORE is free to select locations regardless of infrastructure feasibility or operational accessibility. The constrained SCORE variants reveal the practical cost of enforcing infrastructure feasibility. For Capella, both SCORE (Lat-Const) and SCORE (Infra-Const) perform slightly below unconstrained SCORE but remain close to the Teleports IP baseline, as Capella’s mid-latitude orbits are well-served by infrastructure-accessible locations. Notably, SCORE (Infra-Const) tracks closely with the Teleports line, suggesting that for mid-latitude constellations, restricting placement to infrastructure-accessible regions yields performance comparable to existing teleport networks. For ICEYE, the latitude constraint incurs a more meaningful performance penalty. At small network sizes, constrained SCORE variants perform comparably to fixed-site methods, as the most geometrically valuable polar locations are excluded from the search space. For larger networks, however, SCORE (Lat-Const) and SCORE (Infra-Const) consistently outperform both Teleports and KSAT, achieving gains of 2–5% over Teleports and 2–5% over KSAT at 20 stations. This compares to unconstrained SCORE’s 8–15% improvement, indicating that infrastructure constraints reduce but do not fully eliminate the performance advantage of free placement. These results highlight a key tradeoff for operators: infrastructure constraints reduce the performance ceiling of free placement, particularly for high-inclination constellations, but provide more physically attainable site placements that still can outperform fixed-site selections at larger network sizes. 2. Coordinate Location Comparisons We now perform a detailed evaluation of optimized ground station coordinates for a 10-station network of both the Capella Space (Figure 7) and ICEYE (Figure 8) constellations. In both figures, optimized ground station locations for unconstrained SCORE, KSAT IP solutions, and Teleport IP solutions are displayed on a global map. The figures display optimized station locations on global maps, with coordinates listed in Table 5 and Table 6. Translucent circles represent the communication coverage cones for each ground station, calculated for a satellite altitude of 525 km and a minimum elevation angle of 10◦ . The maps indicate that orbital inclination strongly influences ground station network design. For the ICEYE constellation, which follows polar or sun-synchronous orbits, optimal ground stations cluster at mid-to-high latitudes to maximize contact opportunities. In contrast, the Capella Space constellation, with mid-inclination orbits, has stations

20

Total Data Downlinked (PB) over Topt

Total Data Downlinked (PB) over Topt

70 60 50 40 30

SCORE SCORE (Lat-Const)

20

SCORE (Infra-Const) Teleports KSAT

10

250

200

150

100

SCORE SCORE (Lat-Const)

50

SCORE (Infra-Const) Teleports KSAT

1 2 3 4 5 7 10 15 20 Ground Stations Selected for CAPELLA Space Constellation

1 2 3 4 5 7 10 15 20 Ground Stations Selected for ICEYE Constellation

(a) Full Capella Space Constellation

60 58 56 54

(b) Full ICEYE Constellation

Zoomed View at GS = 15 for CAPELLA Space

Zoomed View at GS = 15 for ICEYE

SCORE SCORE (Lat-Const)

220

SCORE SCORE (Lat-Const)

SCORE (Infra-Const) Teleports KSAT

210

SCORE (Infra-Const) Teleports KSAT

52

200

50 48

190

15

(c) Zoomed View of Capella Constellation

15

(d) Zoomed View of ICEYE Constellation

Fig. 6 Comparison of full constellation views (top) and their corresponding zoomed-in regional views (bottom) for Capella and ICEYE constellations, illustrating free-placement vs. fixed-site ground station selection.

90 60 30 0 −30

0 18

0 15

0 12

90

60

30

0

0 −3

0 −6

0 −9

20 −1

50

−1

80

−90

SCORE KSAT Teleports

−1

−60

Fig. 7 Optimization over Capella Space Constellation showing final ground station locations for a network of 10 stations. Solutions include KSAT stations (red), Teleport sites (purple), and SCORE free-placement (blue).

21

90 60 30 0 −30

0 18

0 15

0 12

90

60

30

0

0 −3

0 −6

0 −9

−1

50 −1

−1

80

−90

20

SCORE KSAT Teleports

−60

Fig. 8 Optimization over ICEYE Constellation showing final ground station locations for a network of 10 stations. Solutions include KSAT stations (red), Teleport sites (purple), and SCORE free-placement (blue).

clustered mainly in mid-latitude regions. While latitude primarily determined how often a ground station can see satellites, longitude was found to affect the timing of satellite passes. A wide spread in longitude seemed to help provide more continuous and evenly distributed downlink opportunities over time. Overall, the combination of latitudinal clustering and longitudinal distribution in the optimized networks demonstrates how both free-placement and fixed-site strategies balance maximizing pass frequency with maintaining temporal coverage. Table 5 summarizes metrics for Capella Space ground stations, showing SCORE downlinked 42.43 PB over 𝑇𝑜 𝑝𝑡 , outperforming KSAT’s 37.44 PB and WTA Teleports’ 41.24 PB. For ICEYE (Table 6), SCORE achieves 165.20 PB, exceeding KSAT’s 151.06 PB and Teleports’ 148.25 PB by 9–11%. SCORE consistently selects more diverse, strategically positioned sites in high-latitude or remote areas, leading to substantial data downlink improvements. These results suggest that expanding ground networks beyond traditional infrastructure hubs to underutilized or remote locations can yield meaningful performance gains. However, practical deployment considerations, including backhaul availability and power infrastructure, remain important factors in translating these geometric gains into operational improvements. Table 5 Comparison of KSAT, Teleport, and SCORE ground stations for Capella: 10 stations placed, with SCORE achieving highest data downlink in 𝑇opt . KSAT Lon.

Lat.

−118.15 −84.26 −25.13 22.69 −70.87 27.71 143.45 −156.45 168.38 115.34

33.82 32.95 36.99 38.82 −52.94 −25.89 42.60 20.82 −46.53 −29.01

Location Long Beach, USA Thomaston, USA Azores, Portugal Thermopylae, Greece Punta Arenas, Chile Hartebeesthoek, S. Africa Hokkaido, Japan Maui, USA Awarua, New Zealand Mingenew, Australia

37.44 PB Downlinked

Teleports Lon.

Lat.

−111.95 −52.78 −3.79 44.90 −58.31 18.72 141.26 −157.87 168.38 115.94

40.78 47.56 40.40 41.70 −34.69 −34.03 43.17 21.31 −46.53 −31.88

Location Salt Lake City, USA St. John’s, Canada Madrid, Spain Tbilisi, Georgia Buenos Aires, Argentina Cape Town, South Africa Hokkaido, Japan Honolulu, USA Awarua, New Zealand Perth, Australia

41.24 PB Downlinked

22

SCORE Lon.

Lat.

−105.98 −59.88 15.10 60.97 −62.77 51.73 107.87 −177.38 148.34 1.92

41.52 43.93 41.64 40.84 −40.30 −46.28 41.04 28.21 −40.39 −89.94

Location Cheyenne, USA Nova Scotia, Canada Foggia, Italy Arap, Turkmenistan Bahía Blanca, Argentina Alfred Faure, TAAF Bayannur, China Sand Island, USA Tasmania, Australia South Pole

43.43 PB Downlinked

3. Site Feasibility and Infrastructure Proximity Analysis To directly address the practical feasibility of free-placement solutions, we analyze the distance between optimized ground station locations and the nearest existing teleport infrastructure in Figure 2b for each SCORE variant. Tables 7 and 8 report coordinates and infrastructure distances for all three variants across both constellations. Unconstrained SCORE selects geometrically optimal sites averaging 1,011 km (Capella) and 860 km (ICEYE) from the nearest teleport infrastructure, including remote locations such as the South Pole (−89.94◦ N) and high-Arctic sites that lack the fiber backhaul and power infrastructure required for commercial EO downlink operations. These results confirm that unconstrained free placement represents a strong empirical performance benchmark that can quantify the substantial gains achievable through unconstrained network placement. The latitude- and infrastructure-constrained variants reveal how much of this gain is recoverable under realistic deployment assumptions. Infrastructure-constrained SCORE reduces mean distance to nearest teleport to 113 km (Capella) and 116 km (ICEYE) — an order-of-magnitude reduction — while retaining 96.5% (41.90/43.43 PB) and 92.0% (152.03/165.20 PB) of unconstrained performance respectively. The latitude-constrained variant retains 98.0% (42.58/43.43 PB) for Capella and 93.5% (154.51/165.20 PB) for ICEYE, at a somewhat higher infrastructure cost (785 km and 810 km mean distance). This demonstrates that substantial geometric gains remain achievable through strategic expansion near existing infrastructure hubs, without requiring deployment at remote or logistically infeasible sites. Thus, free placement defines a strong performance benchmark and identifies strategic regions where targeted infrastructure investment could yield significant gains. The constrained variants provide actionable guidance for network planners operating within real-world deployment constraints, bridging the gap between geometric optimality and operational feasibility. 4. Additional Antennas versus Ground Stations In all evaluations presented in Section V.B, each ground station was constrained to one active satellite link at a time, representing a single-antenna configuration. We further investigated the throughput improvements enabled by equipping stations with multiple antennas. Figures 9 and 10 display 24-hour contact windows for three example ground stations for the Capella Space and ICEYE constellations, highlighting both utilization patterns and the fundamental limits of single-antenna architectures. Both constellations exhibit overlapping contact periods, indicated by red regions. For Capella Space (Figure 9), overlaps recur cyclically across stations. In the ICEYE constellation (Figure 10), overlaps are even more pronounced, especially at high-latitude stations. These overlaps reveal a key limitation of single-antenna stations: simultaneous satellite contacts force prioritization, causing lost downlink opportunities. This challenge intensifies with constellation size, as reflected by substantial red areas in the figures. The analysis suggests that expanding geographic coverage alone may not significantly increase capacity. Instead, equipping existing stations with multiple antennas at strategic locations offers greater throughput gains. Table 9 quantifies this effect: antenna constraints reduce ICEYE throughput by 68% at network size 1 (from 85.29 PB to Table 6 Comparison of KSAT, Teleport, and SCORE ground stations for ICEYE: 10 stations placed, with SCORE achieving the highest data downlink in 𝑇opt . KSAT Lon.

Lat.

−133.72 −148.49 15.63 −8.62 31.11 143.45 −51.72 2.53 −67.52 168.38

68.36 70.26 78.25 71.00 70.37 42.60 64.18 −72.01 −55.05 −46.53

Location Inuvik, Canada Prudhoe Bay, USA Svalbard, Norway Jan Mayen, Norway Vardo, Norway Hokkaido, Japan Nuuk, Greenland Troll, Antarctica Tolhuin, Argentina Awarua, New Zealand

151.06 PB Downlinked

Teleports Lon.

Lat.

−133.55 −148.94 15.39 −21.69 93.51 30.50 −80.33 2.53 −67.12 168.38

68.32 70.32 78.23 64.14 56.25 50.47 41.43 −72.01 −54.51 −46.53

Location Inuvik, Canada Prudhoe Bay, USA Svalbard, Norway Reykjavík, Iceland Zheleznogorsk, Russia Kiev, Russia Greenville PA Troll, Antarctica Tolhuin, Argentina Awarua, New Zealand

148.25 PB Downlinked

23

SCORE Lon.

Lat.

−127.32 175.48 −32.32 −6.74 46.22 102.64 −74.13 66.52 −67.63 158.96

57.78 67.33 83.55 58.13 57.47 66.60 59.81 −90.00 −56.21 −55.14

Location British Columbia, Canada Chukotka, Russia Peary Land, Greenland Scotland, UK Kazan, Russia Chirinda, Russia Tuttusivik, Canada South Pole Tolhuin, Argentina Macquarie Island, Australia

165.20 PB Downlinked

Table 7 Comparison of SCORE, SCORE (Lat-Const), and SCORE (Infra-Const) ground station coordinates for Capella Space: 10 stations placed. Distance to the nearest teleport sites in Figure 2b is reported for each station. SCORE

SCORE (Lat-Const)

Lon.

Lat.

Dist. (km)

Lon.

Lat.

Dist. (km)

107.87 −105.98 148.34 −177.38 15.10 51.73 1.92 −62.77 −59.88 60.97

41.04 41.52 −40.39 28.21 41.64 −46.28 −89.94 −40.30 43.93 40.84

713 112 292 2063 44 3055 1993 487 383 968

−8.14 68.77 −176.95 173.02 109.72 −109.69 136.68 −62.79 20.00 60.97

41.03 −48.50 51.49 −40.77 41.01 40.50 −36.17 −40.32 −34.91 40.84

Mean dist. (km) Max dist. (km)

1011 3055

Mean dist. (km) Max dist. (km)

43.43 PB Downlinked

SCORE (Infra-Const) Lon.

Lat.

Dist. (km)

259 4155 750 176 558 118 231 490 144 968

−159.61 −166.29 140.89 −58.88 136.68 −124.39 19.53 −81.38 −4.41 44.52

21.34 53.35 41.99 −33.72 −36.17 40.25 −34.77 40.98 51.00 41.45

143 61 106 78 231 239 102 11 119 42

785 4155

Mean dist. (km) Max dist. (km)

113 239

42.58 PB Downlinked

41.90 PB Downlinked

27.55 PB) and 53% at network size 20 (from 553.81 PB to 260.89 PB). Capella Space is less impacted, with 7% and 12% penalties at network sizes 1 and 20 respectively, reflecting orbital and contact pattern differences. These results demonstrate that for polar-orbit constellations like ICEYE, multi-antenna stations can more than double throughput without increasing geographic footprint, a critical design consideration for satellite network infrastructure.

Contact No Contact (Gap) Overlapping Contact

Capella 5 - GS 3 Capella 4 - GS 3 Capella 3 - GS 3 Capella 2 - GS 3 Capella 1 - GS 3 Capella 5 - GS 2 Capella 4 - GS 2 Capella 3 - GS 2 Capella 2 - GS 2 Capella 1 - GS 2 Capella 5 - GS 1 Capella 4 - GS 1 Capella 3 - GS 1 Capella 2 - GS 1 Capella 1 - GS 1 04−01 18

04−01 21

04−02 00

04−02 03 04−02 06 04−02 09 Time in MM-DD-HH

04−02 12

04−02 15

04−02 18

Fig. 9 Contact windows for three optimized SCORE ground stations on Capella: blue (contact), gray (no contact), red (overlapping contacts from multiple satellites).

VI. Conclusion In this work, we introduced SCORE, an algorithm for ground station network design that supports free-placement location optimization. SCORE employs a two-step process, sequential coordinate selection followed by sequential cyclic 24

Table 8 Comparison of SCORE, SCORE (Lat-Const), and SCORE (Infra-Const) ground station coordinates for ICEYE: 10 stations placed. Distance to the nearest teleport sites in Figure 2b is reported for each station. SCORE

SCORE (Lat-Const)

Lon.

Lat.

Dist. (km)

Lon.

Lat.

Dist. (km)

−127.32 66.52 −32.32 −74.13 158.96 102.64 −8.13 175.48 46.22 −67.63

57.78 −90.00 85.35 59.81 −55.14 66.60 58.19 67.33 57.47 −56.21

379 2000 1033 1110 1161 953 357 917 500 192

2.73 11.69 −113.62 166.56 150.00 82.46 −67.20 −7.35 70.11 −164.28

−72.5 79.07 70.74 −48.26 72.27 61.40 −56.25 52.96 −49.97 53.45

Mean dist. (km) Max dist. (km)

860 2000

Mean dist. (km) Max dist. (km)

165.20 PB Downlinked

SCORE (Infra-Const) Lon.

Lat.

Dist. (km)

58 123 816 236 1794 637 193 90 4057 156

2.53 166.54 12.69 148.49 −3.82 26.11 −133.58 139.11 −61.63 −67.20

−72.01 −48.27 79.07 70.26 52.06 70.20 69.02 36.61 53.08 −56.25

0 238 110 15 67 270 78 102 90 193

810 4057

Mean dist. (km) Max dist. (km)

116 270

154.51 PB Downlinked

152.03 PB Downlinked

Contact No Contact (Gap) Overlapping Contact

ICEYE 5 - GS 3 ICEYE 4 - GS 3 ICEYE 3 - GS 3 ICEYE 2 - GS 3 ICEYE 1 - GS 3 ICEYE 5 - GS 2 ICEYE 4 - GS 2 ICEYE 3 - GS 2 ICEYE 2 - GS 2 ICEYE 1 - GS 2 ICEYE 5 - GS 1 ICEYE 4 - GS 1 ICEYE 3 - GS 1 ICEYE 2 - GS 1 ICEYE 1 - GS 1 04−01 18

04−01 21

04−02 00

04−02 03 04−02 06 04−02 09 Time in MM-DD-HH

04−02 12

04−02 15

04−02 18

Fig. 10 Contact windows for three optimized SCORE ground stations on the first five ICEYE satellites: blue (contact), gray (no contact), red (overlapping contacts from multiple satellites).

refinement, to effectively address high-dimensional optimization challenges, enabling the design of larger ground station networks to support increasingly complex satellite constellations. We conducted a comprehensive evaluation of SCORE against both single-step free-placement optimization methods, such as differential evolution, and fixed-site selection optimization IP baselines using real-world infrastructure from KSAT and the World Teleport Association. In the free-placement setting, SCORE consistently outperformed DE in both computational efficiency and solution quality. Across our studies, SCORE required 5× fewer function evaluations than DE to reach convergence, while achieving 13% higher data downlink throughput. These results highlight that SCORE’s sequential coordinate selection and local refinement strategy effectively overcomes DE’s limitations in fine-grained exploitation, as our method not only converges faster but also finds higher-quality station placements. Our results

25

Table 9 Comparison of total data downlinked between single-antenna and multi-antenna SCORE for supporting Capella Space and ICEYE constellations. SCORE Capella Space

ICEYE

Network Size

No Antenna Constraint

Antenna Constraint

No Antenna Constraint

Antenna Constraint

(1) (2) (3) (4) (5) (7) (10) (15) (20)

5.04 PB 10.08 PB 15.08 PB 19.94 PB 25.10 PB 34.08 PB 46.35 PB 66.37 PB 82.15 PB

4.70 PB 9.23 PB 13.68 PB 18.04 PB 22.76 PB 30.78 PB 42.41 PB 58.55 PB 72.44 PB

85.29 PB 161.76 PB 200.63 PB 239.48 PB 259.67 PB 310.95 PB 382.38 PB 470.65 PB 553.81 PB

27.55 PB 53.75 PB 70.70 PB 87.75 PB 99.65 PB 125.99 PB 165.20 PB 214.42 PB 260.89 PB

also underscore that the overall composition of station locations plays a more critical role than any individual site in achieving robust satellite communication and data downlink performance. When compared to fixed-site selection optimizations, where IP-based formulations yield optimal solutions within a limited set of locations, SCORE’s unconstrained global search establishes a strong empirical performance benchmark, achieving up to 15% greater total data downlink under flexible siting conditions. This quantifies the geometric performance potential beyond pre-existing site constraints and motivates strategic infrastructure expansion in high-value regions. Infrastructure-constrained SCORE retains over 92% of this gain while remaining near existing fiber and power infrastructure, demonstrating that substantial improvements are achievable within realistic deployment constraints. We also examined the trade-offs between adding antennas to existing stations and expanding the network with strategically placed new ground stations. SCORE offers a practical solution for satellite operators aiming to optimize their ground infrastructure, while also opening new directions for research in scalable and adaptive network design. While developed for free-placement ground station optimization, SCORE’s approach is broadly applicable to other problems requiring structured search over complex configuration spaces. As satellite constellations expand in both size and complexity, methods like SCORE will be essential for efficiently evaluating and designing ground station networks that can effectively support increasing data demands and evolving mission requirements of growing satellite constellations. The complete codebase supporting this study is publicly available at: https://github.com/sisl/loc-gsopt.

Appendix A. Alternative Objective Function: Mean Gap-Time Minimization Another possible objective includes mean gap time minimization, which seeks to reduce the average time between consecutive contacts in order to prevent prolonged communication gaps for any satellite. For an EO mission, this encourages a more consistent communication cadence across the constellation by reducing the typical downtime between uplink and downlink opportunities. For the purpose of computing mean gap times, only feasible contact opportunities are considered (i.e., for each contact pair (𝑐, 𝐶), we set 𝑐 = 1 to indicate that every feasible contact is included in the gap calculation). This ensures the gap metric measures the inherent communication cadence of the network, independent of scheduling decisions. Let C𝑆 ⊆ CL, S be the contacts for a satellite 𝑆 ∈ S supported by ground network L. We assume all contacts (𝑐, 𝐶) ∈ C𝑆 are ordered in ascending order by the start time 𝐶 start . This ordering allows us to define start − 𝐶 end , where 𝐶 end ≤ 𝐶 start . If 𝐶 end ≥ 𝐶 start , the previous contact window communication gaps as the difference 𝐶𝑖+1 𝑖 𝑖 𝑖 𝑖+1 𝑖+1 overlaps with or ends exactly at the start of the next contact window, resulting in continuous communication with no gaps between contacts. For the length of all contacts 𝑖 = 1, . . . , |C𝑆 | − 1, we define the positive gaps lengths 𝑔gaps (L, 𝑆) for each satellite as: start (19) 𝑔gaps (L, 𝑆) = max(0, 𝐶𝑖+1 − 𝐶𝑖end )

26

from contacts generated with a given candidate network L. The number of real gaps for satellite 𝑆 is then defined as ∑︁ 𝑛gaps (L, 𝑆) = 1{𝑔gaps (L, 𝑆) > 0} (20) 𝑆∈S

where only gaps larger than zero seconds are counted. Using both Equation (19) and Equation (20), the per-satellite mean gap time is then calculated as Í 𝑔gaps (L, 𝑆)     𝑆∈S , 𝑛gaps (L, 𝑆) > 0 𝑛gaps (L, 𝑆) 𝑔¯ gaps (L, 𝑆) = (21)   0 , 𝑛gaps (L, 𝑆) = 0  for both when 𝑛 is non-zero, or when there exist no gaps within the mission window. This ensures the mean is computed over actual communication gaps, excluding overlapping contacts. The full objective function then becomes 𝑓gap (L, S) =

𝑇opt 1 ∑︁ 𝑔¯ gaps (L, 𝑆) 𝑇sim |S|

(22)

𝑆∈S

which represents the mean length of communication gap times for each satellite, averaged over the entire constellation. 𝑇opt Once again, we weight the objective by 𝑇sim to appropriately reflect the expected total downtime over the entire mission. B. Synthetic Walker-Star Constellations We outline how to calculate several parameters necessary to define a synthetic Walker-Star constellation for simulation purposes, including the right ascension of the ascending node for each orbital plane, the mean anomaly of each satellite to introduce phasing, and a unique indexing scheme to identify satellites within and across planes. The Right Ascension of the Ascending Node (Ω) for each plane 𝜌 ∈ {1, . . . , 𝑁𝜌 } is uniformly distributed across 360° as Ω𝜌 = 𝜌 ·

360 . 𝑁𝜌

(23)

Each satellite 𝑆 ∈ S in this Walker-Star constellation is uniquely indexed by its orbital plane 𝜌 and its in-plane index 𝑠, written as 𝑆 𝜌,𝑠 . To introduce phasing between planes, we define the mean anomaly 𝑀 of satellite 𝑆 𝜌,𝑠 as 𝑀𝜌,𝑠 =

360 720 𝜌 𝑠+ , 𝑠per plane |S|

(24)

where 𝑠 = 1, . . . , 𝑠per plane is the in-plane satellite index and 𝑠per plane is the number of satellites per orbital plane. The first term in Equation (24) governs the relative placement of satellites within each plane, while the second term introduces inter-plane phasing. For this synthetic Walker-Star constellation, since 𝑠per plane = 1, the in-plane spacing term vanishes, and 𝑀𝜌,𝑠 is determined entirely by inter-plane phasing. C. Surrogate Optimization Analysis As mentioned in section II, long-term orbital propagation cannot reliably predict contact opportunities with sufficient accuracy over multi-year timespans, as stochastic perturbations experienced by LEO spacecraft, specifically solar radiation pressure and atmospheric drag, accumulate over time. To capture the cyclic nature of satellite contacts across orbital cycles, we follow prior work [29, 45] and adopt a 7–10 day surrogate window as a statistically representative approximation of long-term contact distributions. To validate this assumption, we performed a parameter sweep of 5,250 configurations across varying altitudes, inclinations, and ground station locations, detailed in Table 10. For each configuration, contacts were simulated at 365 discrete durations (1-day, 2-day, . . . , 365-day windows), yielding a running mean of contacts/day and duration/day at each duration relative to the 365-day reference. Of the 5,250 configurations, 3,644 yielded non-zero contacts and are included in the analysis. Figure 11 presents running means of contacts/day and duration/day normalized to 365-day references, shown over the first 100 days for clarity. The running mean contact statistics converge steadily toward the full-year mean as simulation

27

Table 10

Surrogate validation parameter sweep.

Parameter

Range

Altitude Inclination KSAT locations Simulation durations

300–1000 km (50 km steps) 0–90◦ (10◦ steps) Stations from the 35 ground station set 1–365 days (1-day steps)

Total configs

5,250 (3,644 non-zero contacts)

duration increases. A 7-day surrogate window estimates annual mean contact frequency and duration within 5% for over 70% of configurations tested, rising to over 95% with a longer window of 20–25 days. Given that propagation accuracy degrades at longer horizons due to accumulated stochastic perturbations, we adopt 7–10 days as a practical surrogate duration balancing statistical representativeness with propagation fidelity. These results substantially extend the preliminary analysis of Eddy et al. [29], providing large-scale empirical justification for this assumption.

1.8

5–95th pct Day 7 365-day mean (reference) ±5

1.6 1.4

Running mean / 365-day mean

Running mean / 365-day mean

1.8

Running Mean Contacts/Day (normalized to 365-day mean)

1.2 1.0 0.8 0.6

Running Mean Contact Duration/Day (normalized to 365-day mean)

% of Configurations Whose Running Mean Is Within 5% of 365-day Mean 100

5–95th pct Day 7 365-day mean (reference) ±5

1.6 1.4

80 60

1.2 1.0

40

0.8

Contacts Duration Day 7 95

20

0.6 20

40

60

Simulation Duration (days)

80

100

20

40

60

Simulation Duration (days)

80

100

0

20

40

60

80

100

Simulation Duration (days)

Fig. 11 Convergence of running mean contact statistics to the 365-day reference across 3,644 configurations. The red dashed line marks the adopted 7-day surrogate window.

D. Differential Evolution Hyperparameter Sweep To ensure a fair comparison between SCORE and Differential Evolution (DE), we conducted a hyperparameter sweep over 9 configurations for our largest synthetic Walker-Star constellation experiment in Section V.A. The sweep spanned crossover rate 𝐶 𝑅 ∈ {0.3, 0.5, 0.9}, mutation scale 𝐹 ∈ {0.5, 0.6, 0.7, 0.9}, population size 𝑁 𝑃 ∈ {5, 10, 15}, and mutation strategy ∈ {best1bin, rand1bin}. We note the final selected configurations in Table 11. These ranges were selected to cover both conservative and aggressive exploration strategies, following standard recommendations from the DE literature [52]. We note that while the majority of DE configurations can eventually approach SCORE’s solution quality given sufficient function evaluations, identifying well-performing hyperparameters itself requires additional tuning effort. The convergence behavior of all configurations is shown in Figure 12. Despite this broad sweep, no configuration matched SCORE’s convergence speed, with most requiring 4–5× more function evaluations to reach comparable solution quality, and two configurations (Runs 1 and 5) failing to converge to competitive solutions entirely. This hyperparameter sensitivity becomes increasingly problematic as ground station network sizes grow, since each function evaluation requires propagating satellite orbits and computing access windows over the full simulation window, making exhaustive parameter sweeps computationally prohibitive at scale. In contrast, SCORE requires no hyperparameter tuning beyond the choice of local optimizer, and consistently converges to near-optimal solutions with modest growth in function evaluations as network size increases.

28

3.5 3.0

Run 1: NP=5d, F=0.5, CR=0.3, best1bin

Run 2: NP=10d, F=0.7, CR=0.9, rand1bin

Run 3: NP=10d, F=0.9, CR=0.3, best1bin

Run 4: NP=15d, F=0.7, CR=0.3, rand1bin

Run 5: NP=5d, F=0.9, CR=0.9, rand1bin

Run 6: NP=10d, F=0.5, CR=0.9, rand1bin

Run 7: NP=15d, F=0.9, CR=0.9, rand1bin

Run 8: NP=10d, F=0.7, CR=0.5, best1bin

Run 9: NP=10d, F=0.6, CR=0.6, best1bin

2.0 1.0

Data Downlinked in Topt (PB)

0.0

3.5 3.0 2.0 1.0 0.0

3.5 3.0

DE Run SCORE Avg

2.0 1.0 0.0 0

5000 10000 15000

0

5000 10000 15000

0

5000 10000 15000

Number of Function Evaluations Fig. 12 Convergence comparison of SCORE versus 9 DE hyperparameter configurations: SCORE converges 4 to 5× faster, with several DE runs failing to converge within the evaluation budget.

29

Table 11 Differential Evolution hyperparameter configurations tested in parameter sweep. Population size is scaled by problem dimensionality 𝑑 (e.g., 5𝑑 means 5 × 𝑑 where 𝑑 = 2𝑛 for 𝑛 ground stations). Run

Pop. Size

𝐹

CR

Strategy

Purpose

1 2 3 4 5 6 7 8 9

5𝑑 10𝑑 10𝑑 15𝑑 5𝑑 10𝑑 15𝑑 10𝑑 10𝑑

0.5 0.7 0.9 0.7 0.9 0.5 0.9 0.7 0.6

0.3 0.9 0.3 0.3 0.9 0.9 0.9 0.5 0.6

best1bin rand1bin best1bin rand1bin rand1bin rand1bin rand1bin best1bin best1bin

Fast convergence, local exploitation Balanced exploration High mutation, low mixing Conservative large population Chaotic, high exploration Baseline (from original experiments) Maximum exploration, large pop Moderate exploration Moderate balanced configuration

E. SCORE Cyclic Refinement Ablation Study We provide an ablation study of the cyclic refinement step in the SCORE algorithm, addressing two questions: (1) whether the cyclic refinement step provides meaningful performance gains over greedy-only placement, and (2) whether randomizing the ground station selection order within the refinement step affects solution quality. 1. Effect of Cyclic Refinement The SCORE algorithm consists of two phases: an initial greedy placement phase, in which ground stations are selected sequentially to maximize marginal coverage, followed by a cyclic refinement phase, in which each station is re-optimized while holding all others fixed, repeated for 𝐸 max cycles. To assess the contribution of the refinement phase, we compare convergence behavior with and without cyclic refinement enabled across the four-satellite, four-station Walker-Star scenario described in Section IV.A.1. As shown in Figure 13, the vertical dashed line at approximately 800 function evaluations marks the boundary between the greedy placement phase and the onset of cyclic refinement. The greedy phase alone achieves a substantial fraction of the final solution quality; however, the cyclic refinement phase consistently yields additional improvement, particularly in runs where the greedy placement produced a suboptimal initial configuration. This demonstrates that the second phase is not redundant and contributes measurably to final solution quality. 2. Effect of Randomized Selection Order We additionally investigate whether randomizing the order in which ground stations are selected for re-optimization during the cyclic refinement phase affects performance. For each of nine independent trials with distinct random seeds, we run SCORE with both a fixed deterministic ordering (𝑖 = 0, 1, . . . , 𝑛 − 1) and a randomly shuffled ordering, where the shuffle is re-drawn at each refinement cycle using a seeded random number generator. Figure 13 shows the convergence curves for both conditions across all nine trials. In the majority of runs, the randomized and deterministic orderings converge to comparable final solutions, with neither consistently outperforming the other. Some runs exhibit transient differences in convergence trajectory, but these differences do not persist to the final solution. This suggests that SCORE’s performance is robust to the specific selection order used in the refinement step, and that the deterministic ordering used in the primary results is not a source of systematic bias.

Funding Sources This research was supported by the Hertz Foundation and the National Defense Science and Engineering Graduate (NDSEG) Fellowship Program.

30

Run 1:

Data Downlinked in Topt (PB)

4 3 2 1 0

Run 2:

Refinement begins

Refinement begins

Run 4:

4 3 2 1 0

Refinement begins

Run 8:

Refinement begins

1500

Run 6:

Refinement begins

Run 7:

0

Refinement begins

Run 5:

Refinement begins

4 3 2 1 0

Run 3:

Run 9:

SCORE Random SCORE No Random

Refinement begins

3000

4500

0

1500

3000

4500

0

1500

3000

4500

Number of Function Evaluations Fig. 13 SCORE solution quality is robust to cyclic refinement ordering: randomized and deterministic variants converge comparably across nine independent Walker-Star trials.

References [1] Eftimiades, N., Small Satellites: The Implications for National Security, Atlantic Council, Scowcroft Center for Strategy and Security, 2022. [2] Falle, A., Wright, E., Boley, A., and Byers, M., “One Million (Paper) Satellites,” Science, Vol. 382, No. 6667, 2023, pp. 150–152. https://doi.org/10.1126/science.adi4639. [3] Karacalioglu, A. G., and Stupl, J., “The Impact of New Trends in Satellite Launches on the Orbital Debris Environment,” International Association for the Advancement of Space Safety Conference Safety First, Safety for All, Melbourne, FL, United States, 2016. URL https://ntrs.nasa.gov/citations/20160011184. [4] Chen, Y., Ma, X., and Wu, C., “The Concept, Technical Architecture, Applications and Impacts of Satellite Internet: A Systematic Literature Review,” Heliyon, Vol. 10, No. 13, 2024. https://doi.org/10.1016/j.heliyon.2024.e33793. [5] Rasila, T., and Ojala, A., “National Regulation of Satellite Ground Stations: A Global Comparison,” Space Business: Emerging Theory and Practice, edited by A. Ojala and W. W. Baber, Springer Nature, Singapore, 2024, pp. 195–217. https://doi.org/10.1007/978-981-97-3430-6_8. [6] Wilkinson, R., Mleczko, M., Brewin, R., Gaston, K., Mueller, M., Shutler, J., Yan, X., and Anderson, K., “Environmental Impacts of Earth Observation Data in the Constellation and Cloud Computing Era,” Science of the Total Environment, Vol. 909, 2024, p. 168584. https://doi.org/10.1016/j.scitotenv.2023.168584.

31

[7] Aragon, B., Houborg, R., Tu, K., Fisher, J. B., and McCabe, M., “CubeSats Enable High Spatiotemporal Retrievals of Crop-Water Use for Precision Agriculture,” Remote Sensing, Vol. 10, No. 12, 2018, p. 1867. https://doi.org/10.3390/rs10121867. [8] Nguyen, T. T., Hoang, T. D., Pham, M. T., Vu, T. T., Nguyen, T. H., Huynh, Q.-T., and Jo, J., “Monitoring Agriculture Areas with Satellite Images and Deep Learning,” Applied Soft Computing, Vol. 95, 2020, p. 106565. https://doi.org/10.1016/j.asoc. 2020.106565. [9] Chuvieco, E., Aguado, I., Salas, J., García, M., Yebra, M., and Oliva, P., “Satellite Remote Sensing Contributions to Wildland Fire Science and Management,” Current Forestry Reports, Vol. 6, No. 2, 2020, pp. 81–96. https://doi.org/10.1007/s40725-02000116-5. [10] Sheffield, J., Wood, E. F., Pan, M., Beck, H., Coccia, G., Serrat-Capdevila, A., and Verbist, K., “Satellite Remote Sensing for Water Resources Management: Potential for Supporting Sustainable Development in Data-Poor Regions,” Water Resources Research, Vol. 54, No. 12, 2018, pp. 9724–9758. https://doi.org/10.1029/2018WR023903. [11] Tao, B., Masood, M., Gupta, I., and Vasisht, D., “Transmitting, Fast and Slow: Scheduling Satellite Traffic Through Space and Time,” International Conference on Mobile Computing and Networking, Association for Computing Machinery, Madrid, Spain, 2023, pp. 1–15. https://doi.org/10.1145/3570361.3592521. [12] Manavalan, “NISAR Real Time Data Processing – A Simple and Futuristic View,” Big Data, Machine Learning, and Applications, edited by R. Patgiri, S. Bandyopadhyay, M. D. Borah, and D. M. Thounaojam, Springer International Publishing, 2020, pp. 95–101. https://doi.org/10.1007/978-3-030-62625-9_9. [13] Vasisht, D., and Chandra, R., “A Distributed and Hybrid Ground Station Network for Low Earth Orbit Satellites,” Association for Computing Machinery Workshop on Hot Topics in Networks, Association for Computing Machinery, Virtual Event, USA, 2020, pp. 190–196. https://doi.org/10.1145/3422604.3425926. [14] Wang, X., Yu, R., Yang, D., and Xue, G., “Infiltrating the Sky: Data Delay and Overflow Attacks in Earth Observation Constellations,” Institute of Electrical and Electronics Engineers International Conference on Network Protocols, 2024, pp. 1–11. https://doi.org/10.1109/ICNP61940.2024.10858559. [15] Linares, L., Vazquez, R., Perea, F., and Galán-Vioque, J., “A Mixed Integer Linear Programming Model for Resolution of the Antenna-Satellite Scheduling Problem,” Institute of Electrical and Electronics Engineers Transactions on Aerospace and Electronic Systems, Vol. 60, No. 1, 2024, pp. 463–473. https://doi.org/10.1109/TAES.2023.3326422. [16] Cho, D.-H., Kim, J.-H., Choi, H.-L., and Ahn, J., “Optimization-Based Scheduling Method for Agile Earth-Observing Satellite Constellation,” Journal of Aerospace Information Systems, Vol. 15, No. 11, 2018, pp. 611–626. https://doi.org/10.2514/1.I010620. [17] Tian, M., Ma, G., Huang, P., Cheng, B., and Li, W., “Optimizing Satellite Ground Station Facilities Scheduling for RSGS: a Novel Model and Algorithm,” International Journal of Digital Earth, Vol. 16, No. 1, 2023, pp. 3949–3972. https://doi.org/10.1080/17538947.2023.2259870. [18] Zilberstein, I., Rao, A., Salis, M., and Chien, S., “Decentralized, Decomposition-Based Observation Scheduling for a Large-Scale Satellite Constellation,” International Conference on Automated Planning and Scheduling, Vol. 34, 2024, pp. 716–724. https://doi.org/10.1609/icaps.v34i1.31535. [19] Fitzgibbon, L., and Shaffer, C., “Simple Scheduling Algorithm for Growing Constellations,” Small Satellite Conference, 2023. URL https://digitalcommons.usu.edu/smallsat/2023/all2023/277. [20] Maule, R., Panigrahy, N. K., Anipeddi, N. L., Dhara, P., Kilbane, D., Hossain, M. Z., Krawec, W. O., Towsley, D., and Wang, B., “Fair and Efficient Scheduling Strategies for Satellite Assisted Quantum Key Distribution Systems,” Institute of Electrical and Electronics Engineers International Conference on Quantum Computing and Engineering, 2024, pp. 1788–1798. https://doi.org/10.1109/QCE60285.2024.00208. [21] Guo, J., Rincon, D., Sallent, S., Yang, L., Chen, X., and Chen, X., “Gateway Placement Optimization in LEO Satellite Networks Based on Traffic Estimation,” Institute of Electrical and Electronics Engineers Transactions on Vehicular Technology, Vol. 70, No. 4, 2021, pp. 3860–3876. https://doi.org/10.1109/TVT.2021.3065994. [22] Liu, S., Wu, T., Hu, Y., Xiao, Y., Wang, D., and Liu, L., “Throughput Evaluation and Ground Station Planning for LEO Satellite Constellation Networks,” International Conference on Space Information Network, Springer, 2019, pp. 3–15. https://doi.org/https://doi.org/10.1007/978-981-15-3442-3_1.

32

[23] Baeza, V. M., Ortiz, F., Lagunas, E., Abdu, T. S., and Chatzinotas, S., “Gateway Station Geographical Planning for Emerging Non-Geostationary Satellites Constellations,” Institute of Electrical and Electronics Engineers Network, Vol. 38, No. 4, 2023, pp. 158–165. https://doi.org/10.1109/MNET.2023.3321531. [24] Baeza, V. M., Ortiz, F., Lagunas, E., Abdu, T. S., and Chatzinotas, S., “Multi-Criteria Ground Segment Dimensioning for NonGeostationary Satellite Constellations,” Joint European Conference on Networks and Communications & 6G Summit, Institute of Electrical and Electronics Engineers, 2023, pp. 252–257. https://doi.org/10.1109/EuCNC/6GSummit58263.2023.10188237. [25] Abe, Y., Ortiz, F., Lagunas, E., Baeza, V. M., Chatzinotas, S., and Tsuji, H., “Optimizing Satellite Network Infrastructure: A Joint Approach to Gateway Placement and Routing,” Institute of Electrical and Electronics Engineers Vehicular Technology Conference, Institute of Electrical and Electronics Engineers, 2024, pp. 1–6. https://doi.org/10.1109/VTC2024-Spring62846.2024.10683455. [26] Chen, Q., Yang, L., Liu, X., Guo, J., Wu, S., and Chen, X., “Multiple Gateway Placement in Large-Scale Constellation Networks with Inter-Satellite Links,” International Journal of Satellite Communications and Networking, Vol. 39, No. 1, 2021, pp. 47–64. https://doi.org/https://doi.org/10.1002/sat.1353. [27] Del Portillo, I., Sanchez, M., Cameron, B., and Crawley, E., “Architecting the ground segment of an optical space communication network,” Institute of Electrical and Electronics Engineers Aerospace Conference, Institute of Electrical and Electronics Engineers, 2016, pp. 1–13. https://doi.org/10.1109/AERO.2016.7500803. [28] Del Portillo, I., Sanchez, M., Cameron, B., and Crawley, E., “Optimal Location of Optical Ground Stations to Serve LEO Spacecraft,” Institute of Electrical and Electronics Engineers Aerospace Conference, Institute of Electrical and Electronics, Big Sky, MT, USA, 2017, pp. 1–16. https://doi.org/10.1109/AERO.2017.7943631. [29] Eddy, D., Ho, M., and Kochenderfer, M. J., “Optimal Ground Station Selection for Low-Earth Orbiting Satellites,” Institute of Electrical and Electronics Engineers Aerospace Conference, 2025, pp. 1–13. https://doi.org/10.1109/AERO63441.2025. 11068558. [30] Fuchs, C., Poulenard, S., Perlot, N., Riedi, J., and Perdigues, J., “Optimization and Throughput Estimation of Optical Ground Networks for LEO-Downlinks, GEO-Feeder Links and GEO-Relays,” Free-Space Laser Communication and Atmospheric Propagation, Vol. 10096, Society For Optics & Photonics, 2017, pp. 298–307. https://doi.org/10.1117/12.2254795. [31] Fuchs, C., and Moll, F., “Ground Station Network Optimization for Space-to-Ground Optical Communication Links,” Journal of Optical Communications and Networking, Vol. 7, No. 12, 2015, pp. 1148–1159. https://doi.org/10.1364/JOCN.7.001148. [32] Cui, H., Chen, X., Guo, M., Jiao, Y., Cao, J., and Qiu, J., “A Distribution Center Location Optimization Model Based on Minimizing Operating Costs Under Uncertain Demand with Logistics Node Capacity Scalability,” Physica A: Statistical Mechanics and its Applications, Vol. 610, 2023, p. 128392. https://doi.org/10.1016/j.physa.2022.128392. [33] Sachan, R., Choi, T. J., and Ahn, C. W., “A Genetic Algorithm with Location Intelligence Method for Energy Optimization in 5G Wireless Networks,” Discrete Dynamics in Nature and Society, Vol. 2016, No. 1, 2016, p. 5348203. https://doi.org/10.1155/ 2016/5348203. [34] Ayati, A., Naji, H. R., Hashemi, M. M., and Saffar, M., “Optimizing Location Allocation in Urban Management: A Brief Review,” International Computer Conference, Computer Society of Iran, 2025, pp. 1–7. https://doi.org/10.1109/CSICC65765. 2025.10967426. [35] Kim, J., Yang, H., and Choe, J., “Robust Optimization of the Locations and Types of Multiple Wells using CNN Based Proxy Models,” Journal of Petroleum Science and Engineering, Vol. 193, 2020, p. 107424. https://doi.org/10.1016/j.petrol.2020. 107424. [36] Wang, L., Yao, Y., Luo, X., Daniel Adenutsi, C., Zhao, G., and Lai, F., “A Critical Review on Intelligent Optimization Algorithms and Surrogate Models for Conventional and Unconventional Reservoir Production Optimization,” Fuel, Vol. 350, 2023, p. 128826. https://doi.org/10.1016/j.fuel.2023.128826. [37] Raisanen, L., and Whitaker, R. M., “Comparison and Evaluation of Multiple Objective Genetic Algorithms for the Antenna Placement Problem,” Mobile Networks and Applications, Vol. 10, No. 1, 2005, pp. 79–88. https://doi.org/10.1023/B: MONE.0000048547.84327.95. [38] Guner, A. R., and Sevkli, M., “A Discrete Particle Swarm Optimization Algorithm for Uncapacitated Facility Location Problem,” Journal of Artificial Evolution and Applications, Vol. 2008, No. 1, 2008, p. 861512. https://doi.org/10.1155/2008/861512. [39] Xu, Q., Zhang, L., and Yu, W., “A Localization Method of Ant Colony Optimization in Nonuniform Space,” Sensors, Vol. 22, No. 19, 2022, p. 7389. https://doi.org/10.3390/s22197389.

33

[40] Crawford, L., Cheng, V., Burns, R., and Liu, S., “Near-Optimal Antenna Placement Using Genetic Search,” Symposium on Multidisciplinary Analysis and Optimization, 2000. https://doi.org/10.2514/6.2000-4914. [41] Liu, J.-L., Chou, C.-W., and Chen, C.-M., “Optimising Mobile Base Station Placement Using an enhanced Multi-Objective Genetic Algorithm,” International Journal of Business Intelligence and Data Mining, Vol. 5, No. 1, 2010, pp. 19–42. https://doi.org/10.1504/IJBIDM.2010.030297. [42] Bagchi, T. P., “Near Optimal Ground Support in Multi-Spacecraft Missions: A GA Model and its Results,” Institute of Electrical and Electronics Engineers Transactions on Aerospace and Electronic Systems, Vol. 45, No. 3, 2009, pp. 950–964. https://doi.org/10.1109/TAES.2009.5259176. [43] Gurobi Optimization, LLC, “Gurobi Optimizer Reference Manual,” https://www.gurobi.com, 2024. [44] Lougee-Heimer, R., “The Common Optimization INterface for Operations Research: Promoting Open-Source Software in the Operations Research Community,” International Business Machines Journal of Research and Development, Vol. 47, No. 1, 2003, pp. 57–66. https://doi.org/10.1147/rd.471.0057. [45] Kim, G. R., Eddy, D., Srinivas, V., and Kochenderfer, M. J., “Scalable Ground Station Selection for Large LEO Constellations,” Institute of Electrical and Electronics Engineers Aerospace Conference, Institute of Electrical and Electronics Engineers, 2026. [46] Nelder, J. A., and Mead, R., “A Simplex Method for Function Minimization,” The Computer Journal, Vol. 7, No. 4, 1965, pp. 308–313. https://doi.org/10.1093/comjnl/7.4.308. [47] Powell, M. J., “An Efficient Method for Finding the Minimum of a Function of Several Variables without Calculating Derivatives,” The Computer Journal, Vol. 7, No. 2, 1964, pp. 155–162. https://doi.org/https://doi.org/10.1093/comjnl/7.2.155. [48] Ijaz, S., Hamayun, M. T., Yan, L., and Mumtaz, M. F., “Fractional Order Modeling and Control of Twin Rotor Aero Dynamical System using Nelder Mead Optimization,” Journal of Electrical Engineering and Technology, Vol. 11, No. 6, 2016, pp. 1863–1871. https://doi.org/10.5370/JEET.2016.11.6.1863. [49] Bezerra, M. A., dos Santos, Q. O., Santos, A. G., Novaes, C. G., Ferreira, S. L. C., and de Souza, V. S., “Simplex Optimization: A Tutorial Approach and Recent Applications in Analytical Chemistry,” Microchemical Journal, Vol. 124, 2016, pp. 45–54. https://doi.org/10.1016/j.microc.2015.07.023. [50] Nobahari, H., Zandavi, S. M., and Mohammadkarimi, H., “Simplex Filter: A Novel Heuristic Filter for Nonlinear Systems State Estimation,” Applied Soft Computing, Vol. 49, 2016, pp. 474–484. https://doi.org/10.1016/j.asoc.2016.08.008. [51] Mohsin, A., Alsmadi, Y., Arshad Uppal, A., and Gulfam, S. M., “A Modified Simplex Based Direct Search Optimization Algorithm for Adaptive Transversal FIR Filters,” Science Progress, Vol. 104, No. 2, 2021, p. 00368504211025409. https://doi.org/10.1177/00368504211025409. [52] Storn, R., and Price, K., “Differential Evolution – A Simple and Efficient Heuristic for Global Optimization over Continuous Spaces,” Vol. 11, No. 4, 1997, pp. 341–359. https://doi.org/10.1023/A:1008202821328. [53] Vallado, D., Crawford, P., Hujsak, R., and Kelso, T. S., “Revisiting Spacetrack Report #3,” AIAA/AAS Astrodynamics Specialist Conference, 2006, p. 6753. https://doi.org/10.2514/6.2006-6753.

34

Record · ID 271754 · SHA-256 2b6626aa3863ec2a
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.