Spacecraft Reachable Domain and Its Applications in Orbital Games: A Review and Future Perspectives

1. Introduction

The rapid proliferation of space activities has fundamentally transformed the orbital environment. With the advent of mega-constellations, on-orbit servicing, active debris removal, and increasingly congested space traffic, the need for systematic analysis of spacecraft maneuvering capabilities has never been greater [1,2]. In many operational scenarios, spacecraft do not operate in isolation but must interact with other objects—whether cooperative targets, non-cooperative adversaries, or uncontrolled debris. These interactions raise fundamental questions about interception feasibility, collision risk, and strategic decision-making in shared orbital environments [3,4,5,6].

The concept of the reachable domain (RD) provides a powerful geometric framework for such analyses. The RD is defined as the set of all terminal states attainable from a given initial state under specified maneuverability constraints, capturing the complete maneuvering capability of a spacecraft in a compact set-theoretic form [7,8,9,10]. Unlike traditional trajectory optimization, which seeks a single optimal path [11,12], the RD enables set-based reasoning about interception feasibility, collision risk, and game-theoretic interactions [13]. This set-theoretic perspective shifts the analysis from point-wise trajectory design to region-based reachability assessment, providing a foundation for subsequent developments in uncertain and adversarial environments.

Building upon this set-theoretic foundation, a variety of methods have been developed over the past decade to compute RDs under different assumptions about thrust and dynamics. Early work focused on single-impulse maneuvers under two-body dynamics. Xue et al. [7] established the foundational classification of single-impulse RDs, distinguishing scenarios based on whether the maneuver point and impulse direction are fixed or free. Wen et al. [8,14] subsequently developed analytical methods for precise boundary determination, introducing the algebraic criterion and extremum envelope approaches. The characteristic ellipse approximation [15] further reduced computational cost, enabling real-time applications. For relative motion scenarios, Wen and Gurfil [16] extended the framework to relative RDs under initial state uncertainties, while Zagaris and Hess [17] provided efficient linearized approximations for short-duration maneuvers. More recently, researchers have extended RD computation to complex dynamical regimes. These include perturbed orbits via Jet Transport techniques [18] and low-thrust propulsion in the circular restricted three-body problem[19]. Polynomial-based methods for arbitrary dynamics have also emerged [20], reflecting a growing emphasis on accuracy under realistic conditions.

Unlike impulsive maneuvers, which assume instantaneous velocity changes, continuous thrust propulsion requires different treatment of the RD problem. Wei and Li [10] compared absolute RDs under continuous finite thrust with impulsive counterparts, revealing significant morphological differences. Wang and Jiang [21,22] derived analytical optimal solutions for low-thrust RDs using Pontryagin’s maximum principle. Bowerfind and Taheri [23] developed rapid approximation methods for low-thrust reachable sets in complex dynamics, while Pang and Wen [24] proposed a grid method for non-convex control constraints. For elliptical reference orbits, Lin and Zhang [25] extended continuous-thrust RD analysis to linear time-varying relative dynamics described by the Tschauner-Hempel equations.

The methods discussed above all belong to the spatial-domain category, focusing on the maximum spatial envelope reachable under fuel or energy constraints without restricting flight time. Time-scale methods answer what positions a spacecraft can reach at a specific time or within a time window. Early work by Venigalla et al. [26] derived Delta-V-based conditions for time-constrained pursuit-evasion. Aguilar-Marsillach et al. [27] proposed reachability-based methods using orbital element differences for passive safety guarantees, while Yang et al. [28] introduced a physics-informed deep learning framework for task assessment under time constraints. Building on these foundations, Zhang et al. [29] developed an exact envelope determination method solving Kepler’s equation numerically. This method systematically distinguishes prescribed-time interception from time-window interception. For onboard applications requiring millisecond-level computation, Yang et al. [30] proposed a fast analytical approximation based on series expansion, while the linearized ellipsoidal approximation [31] offers a computationally lighter alternative for short-duration and small-impulse scenarios.

Building upon these computational methods, RD theory has been widely applied to orbital games. In pursuit-evasion games, researchers have used RDs to assess interception feasibility, compute minimum required Delta-V, and design optimal pursuit strategies [32,33,34]. For single-spacecraft pursuit-evasion scenarios, the fundamental condition for guaranteed interception is that the evader’s RD is contained within that of the pursuer. Safety, in contrast, requires that the two domains be disjoint [29,35]. Time-constrained scenarios have been addressed through fixed-time RD, with Zhang et al. [36] proposing analytical impulse thrust strategies that account for mass variation and perception frequency. This framework extends naturally to multi-agent cooperative games. Li et al. [37] studied cooperative encirclement using spacecraft swarms, showing that three chasers in triangular formation provide optimal planar coverage. Li and Zhang [38] extended the analysis to relative RDs in target-centered frames, deriving "capture regions" for chaser positioning. For cooperative interception, Qing and Zhang [39] developed task allocation methods based on RD coverage, while Wang and Zhang [40] introduced "reachable domain compression" for pursuit under sparse rewards. From the learning perspective, Qiao et al. [41] introduced a RD-guided reinforcement learning method for spacecraft game decision-making.

The dual problem of pursuit is threat avoidance and defense, where the focus shifts from capturing an adversary to avoiding capture or protecting valuable assets. Early work by Zhang et al. [42] incorporated trajectory safety constraints into RD analysis. Zhang et al. [35,43] subsequently developed threat avoidance strategies using predictive RDs as threat envelopes. For passive threats such as debris clouds, Wen and Gurfil [44] modeled debris cloud evolution, while Xu et al. [45] proposed zonotope-based methods for collision prediction. Li and Zhang [46,47] formulated collision avoidance as a reachable domain distance maximization problem. For active threats, Zhang et al. [48] introduced the "escape zone" concept for optimal evasion. Cooperative defense has been addressed by Zhang et al. [49], who proposed a multi-spacecraft cooperative guard strategy. Beyond orbital games, RD theory has found applications in asteroid landing and hopping [50,51,52], ground track reachability [53], reentry flight capability assessment [54], orbit protection [55], and cislunar mission analysis [56,57,58].

Despite significant progress, existing reviews have been either narrowly focused on specific methodological aspects or scattered across disparate application domains. Feng et al. [59] provided a review of RD methods from deterministic to uncertain environments, but their treatment of orbital game applications remains limited. Cao et al. [60] reviewed intelligent orbital game technologies, yet the central role of RDs deserves deeper elaboration. A comprehensive review that bridges the unified modeling framework, solution methods, and orbital game applications of reachable domains is currently lacking. To fill this gap, the present paper provides a systematic review that makes three key contributions: (i) a unified three-dimensional classification (spatial, temporal, and informational attributes) that captures the full scope of RD problems; (ii) a structured overview of solution methods categorized by spatial-domain and time-scale approaches, with comparative guidance for method selection; and (iii) a detailed analysis of orbital game applications across pursuit-evasion, cooperative encirclement, and threat avoidance scenarios, highlighting the dual role of deterministic and predictive RDs.

The remainder of this review is organized as follows. Section 2 establishes a unified mathematical framework for the spacecraft RD, introducing three complementary classification dimensions. Section 3 reviews solution methods for computing spacecraft RDs, organized according to the spatial-temporal classification. Section 4 examines applications of RD theory to orbital games through three representative scenarios: pursuit-evasion, cooperative encirclement, and threat avoidance. Section 5 identifies limitations of existing approaches and proposes future research directions. Finally, Section 6 concludes the review with a summary of key contributions and a forward-looking perspective.

2. Problem Description and Unified Modeling of Spacecraft Reachable Domain

In this section, we establish the fundamental mathematical description of the spacecraft RD and propose a unified modeling framework that accommodates various mission scenarios. The framework is organized around three complementary attributes that capture different facets of the reachability problem. Specifically, we introduce the basic set-theoretic definition of the RD, and then classify it along three dimensions: spatial attributes (absolute and relative), temporal attributes (free-time and fixed-time), and informational attributes (deterministic and predictive). This three-dimensional classification provides the theoretical foundation for the solution methods and orbital game applications discussed in subsequent sections.

2.1. Basic Definitions and Classification

The spacecraft RD is fundamentally defined as the set of all terminal states that can be attained from a given initial state under specified maneuverability constraints. Mathematically, it can be expressed as:

$$\mathcal{R}=\left\{x\in \mathcal{X}|\exists u(·)\in \mathcal{U},\dot{x}=F(x,u,t),x\left(t_0\right)=x_0,c(x,u,t)\le 0\right\}$$

(1)

where x denotes the state vector (typically including position and velocity), $u(·)$ is the control input belonging to admissible control set, $F(·)$ represents the dynamics governing spacecraft motion (e.g., two-body or perturbed orbital dynamics). The constraint function $c(·)$ encapsulates various mission constraints, including path constraints (e.g., keep-out zones), control saturation (e.g., maximum thrust magnitude), and safety boundaries (e.g., minimum approach distance). This unified set-theoretic formulation serves as the basic for all subsequent discussions.

As a geometric description of spacecraft maneuverability, the representation of the RD directly influences both the mathematical formulation and the choice of solution methods. Depending on the analytical dimension and application requirements, RDs can be classified from three complementary perspectives: spatial attributes, temporal attributes, and informational attributes, as illustrated in Figure 1.

2.2. Classification by Spatial Attributes

From the perspective of spatial attributes, the RD can be categorized into absolute RD and relative RD [61,62]. These two categories reflect different physical scenarios and engineering requirements.

The absolute RD refers to the set of all absolute spatial positions that a spacecraft can reach under given initial orbital conditions and maneuverability constraints. It is typically studied within an Earth-centered inertial frame, focusing on the positional reachable range in three-dimensional space. Its mathematical description takes the form:

$${\mathcal{R}}_{\mathrm{abs}}\left(\Delta v_{\max }\right)=\left\{\mathbf{r}\in {\mathbb{R}}^3|\begin{array}{cc}\hfill & \exists f\in [0,2\pi ),\alpha \in [0,2\pi ),\beta \in [-\pi /2,\pi /2],\parallel \Delta \mathbf{v}\parallel \le \Delta v_{\max },\hfill \\ \hfill & \mathrm{s}.\mathrm{t}.\mathbf{r}=\mathbf{g}(f,\alpha ,\beta ;{\mathbf{r}}_0,{\mathbf{v}}_0,\Delta \mathbf{v})\hfill \end{array}\right\}$$

(2)

where ${\mathbf{r}}_0$, ${\mathbf{v}}_0$ denote the initial position and velocity vectors, f is the true anomaly at the maneuver point, $\alpha$ and $\beta$ are impulse direction angles, and $\Delta v_{\max }$ is the maximum available velocity increment. The absolute RD provides an intuitive characterization of a spacecraft’s overall maneuvering capability, making it particularly suitable for capability assessment and target screening in the early stages of mission planning [63,64].

The relative RD describes the set of all possible relative positions that a chasing spacecraft can attain with respect to a target spacecraft (or reference point) under given maneuverability constraints. Following the definition in the literature [16,65], the relative RD can be expressed as:

$${\mathcal{R}}_{rel}\left(t\right)=\left\{\rho \left(t\right)\in {\mathbb{R}}^3|\begin{array}{cc}\hfill & \exists \mathbf{u}(·)\in \mathcal{U},\mathbf{X}\left(t_0\right)\in {\mathcal{X}}_0,\hfill \\ \hfill & \dot{\mathbf{X}}=A\left(t\right)\mathbf{X}+B\left(t\right)\mathbf{u},\rho \left(t\right)=H\mathbf{X}\left(t\right)\hfill \end{array}\right\}$$

(3)

where:

• $\rho \left(t\right)={[x\left(t\right),y\left(t\right),z\left(t\right)]}^T$ is the relative position vector in the target-centered reference frame;
• $\mathbf{X}=[\rho ;\dot{\rho }]\in {\mathbb{R}}^6$ is the augmented relative state vector comprising both relative position and velocity;
• $A\left(t\right)\in {\mathbb{R}}^{6\times 6}$ and $B\left(t\right)\in {\mathbb{R}}^{6\times m}$ are the dynamics and control influence matrices characterizing the relative motion. For circular reference orbits, the Clohessy-Wiltshire (C-W) equations provide a linear time-invariant formulation [66]; for elliptical reference orbits, the Tschauner-Hempel (T-H) equations yield a linear time-varying system;
• ${\mathcal{X}}_0\subseteq {\mathbb{R}}^6$ is the set of admissible initial relative states. In deterministic scenarios, ${\mathcal{X}}_0$ reduces to a singleton $\left\{{\mathbf{X}}_0\right\}$; in the presence of initial state uncertainties or non-cooperative settings, ${\mathcal{X}}_0$ is typically modeled as an ellipsoid or a convex polytope [16];
• $\mathcal{U}$ is the admissible control set, which encapsulates fuel/energy constraints. For impulsive thrust, $\mathcal{U}=\{\Delta \mathbf{v}\mid \parallel \Delta \mathbf{v}\parallel \le \Delta v_{max}\}$; for continuous thrust, $\mathcal{U}=\{{\mathbf{a}}_c(·)\mid \parallel {\mathbf{a}}_c\left(t\right)\parallel \le a_{max}\}$;
• $H=[I_{3\times 3},0_{3\times 3}]\in {\mathbb{R}}^{3\times 6}$ is the position extraction matrix.

This formulation unifies two important perspectives. When the initial state is exactly known (${\mathcal{X}}_0=\left\{{\mathbf{X}}_0\right\}$) and the control inputs are fully specified ($\mathbf{u}(·)$ known), the relative RD reduces to a single trajectory. When only bounds on the control inputs or initial uncertainties are available, the relative RD becomes a set—either a convex envelope (e.g., ellipsoid or zonotope) or a non-convex reachable set—which is precisely the scenario encountered in orbital games where an adversary’s true maneuver cannot be known a priori.

The relative RD has significant applications in close-proximity operations such as rendezvous and docking, formation flying, on-orbit servicing, and non-cooperative target interception [59]. In orbital game scenarios, the relative RD of an adversary serves as a threat envelope that guides the design of safe and effective pursuit, encirclement, or evasion strategies.

2.3. Classification by Temporal Attributes

From the perspective of temporal attributes, RDs are classified into free-time RD and fixed-time RD according to whether the flight time is specified. A further distinction can be made between the RD at a given time (instantaneous) and over a duration (cumulative), as formalized in the time-dependent reachable domain (TDRD) framework [29,31].

The free-time RD represents the set of positions attainable without any restriction on flight time, considering only energy or fuel constraints. This is the most common form in the literature, representing the maximum spatial coverage under given maneuverability constraints [9]. Feng et al. [59] noted that its envelope depends only on the maximum available energy, a property that simplifies the solution process.

The fixed-time RD refers to the set of positions reachable at a specific time instant or within a specified time window. This class is intimately connected with Lambert’s problem, as it must satisfy both position and time constraints. Fixed-time RDs are essential for time-critical missions such as rendezvous at a prescribed epoch or emergency return within a survivable window.

Within the fixed-time framework, Zhang et al. [29] distinguished two subcategories: the RD at a given time $\mathcal{R}\left(t_f\right)$ (exactly at $t_f$) and the RD over a duration $\mathcal{R}\left([t_0,t_f]\right)={\bigcup }_{\tau \in [t_0,t_f]}\mathcal{R}\left(\tau \right)$. This distinction is important in pursuit-evasion scenarios, where the interceptor may have flexibility in interception time [31].

Denoting $\mathcal{R}\left(T\right)$ as the fixed-time RD and ${\mathcal{R}}_{\infty }$ as the free-time RD, the following relationship holds:

$$\mathcal{R}\left(T\right)\subseteq {\mathcal{R}}_{\infty },{\mathcal{R}}_{\infty }=\bigcup _{t_f\ge t_0}\mathcal{R}\left(t_f\right)$$

(4)

This reveals a key insight: the free-time RD is the union of all fixed-time RDs over all possible terminal times, rather than a limiting case at infinite time. In other words, any point in ${\mathcal{R}}_{\infty }$ can be reached within some finite time $t_f$, though that time may be arbitrarily long. As the allowed time increases, the fixed-time RD expands monotonically, and every point in ${\mathcal{R}}_{\infty }$ is eventually covered by $\mathcal{R}\left(T\right)$ for some sufficiently large but finite T.

For small impulsive maneuvers under linearized dynamics, Ma and Zhang [31] proposed an ellipsoidal approximation based on the state transition matrix:

$$\mathcal{R}\left(t_f\right)\approx \left\{\rho \left(t_f\right)=\Phi (t_f,t_0){\rho }_0+\Psi (t_f,t_0)\Delta \mathbf{v}|\parallel \Delta \mathbf{v}\parallel \le \Delta v_{max}\right\}$$

(5)

where $\Phi (t_f,t_0)$ is the state transition matrix that maps the initial relative state to the terminal relative state under zero control, and $\Psi (t_f,t_0)$ is the control influence matrix that maps the impulse vector $\Delta \mathbf{v}$ to the terminal relative position. This linearized approximation is computationally efficient but accurate only for small maneuvers and short propagation times.

Zhang et al. [29] provided a more rigorous formulation without linearization, explicitly enforcing Keplerian propagation:

$$\mathcal{R}\left(t_f\right)=\left\{\mathbf{r}\left(t_f\right)\in {\mathbb{R}}^3|\begin{array}{cc}\hfill & \exists \Delta \mathbf{v},\parallel \Delta \mathbf{v}\parallel \le \Delta v_{max},\mathbf{r}\left(t_f\right)=\mathbf{r}(t_0;{\mathbf{r}}_0,{\mathbf{v}}_0,\Delta \mathbf{v}),\hfill \\ \hfill & \mathrm{Kepler}\text{’}\mathrm{s}\mathrm{equation}\mathrm{satisfied}\mathrm{at}{t}_f\hfill \end{array}\right\}$$

(6)

The distinction between the two classes has profound implications. Free-time RDs are governed by energy-optimality principles, while fixed-time RDs are governed by time-optimality principles. This manifests in solution approaches: free-time problems often reduce to extremum seeking over orbital parameters, whereas fixed-time problems require solving Lambert-type boundary value problems. From an applications perspective, free-time RDs suit early mission design and capability assessment, while fixed-time RDs are indispensable for time-critical operations such as pursuit-evasion and emergency return.

2.4. Classification by Informational Attributes

A distinctive addition to the classification framework, motivated by orbital game scenarios, is the informational attribute perspective. Unlike spatial and temporal attributes, which characterize intrinsic geometric and temporal properties, the informational attribute addresses the epistemic uncertainty associated with control inputs. In non-cooperative settings, the RD assumes two fundamentally different roles: the deterministic RD and the predictive (threat) RD. Figure 2 illustrates this distinction schematically.

The deterministic RD applies to scenarios where the spacecraft’s own control inputs are fully known. Here, the admissible control set $\mathcal{U}$ is precisely characterized (e.g., $\parallel \Delta \mathbf{v}\parallel \le \Delta v_{max}$). The RD serves as a planning tool for trajectory design: given known initial state ${\mathbf{x}}_0$ and control capabilities, the spacecraft can determine exactly which terminal states are achievable. Its computation reduces to solving an optimal control or reachability problem with complete information.

In contrast, the predictive RD applies to scenarios where the target’s control inputs are unknown, and only a bound on its maneuvering capability (e.g., $\Delta v_{max}^a$) is available. This is the typical situation in non-cooperative orbital games: a pursuer knows the evader’s fuel budget but not the exact timing, direction, or magnitude of its maneuvers. The RD then serves as a prediction tool to estimate the set of all possible states that the target could reach [67,68]. Formally:

$${\mathcal{R}}_{\mathrm{pred}}(t_f\mid \Delta v_{max}^a)=\left\{{\mathbf{x}}_a\left(t_f\right)|\begin{array}{cc}\hfill & \exists {\mathbf{u}}_a(·)\in {\mathcal{U}}_a\left(\Delta v_{max}^a\right),{\dot{\mathbf{x}}}_a={\mathbf{f}}_a({\mathbf{x}}_a,{\mathbf{u}}_a,t),\hfill \\ \hfill & {\mathbf{x}}_a\left(t_0\right)={\mathbf{x}}_a^0\hfill \end{array}\right\}$$

(7)

This predictive set defines a threat envelope. For safe operation, the ego spacecraft’s trajectory should be disjoint from this envelope: ${\mathcal{R}}_{\mathrm{ego}}\left(t\right)\cap {\mathcal{R}}_{\mathrm{pred}}\left(t\right)=⌀$. For interception, the pursuer seeks the opposite: ${\mathcal{R}}_{\mathrm{pursuer}}\left(t\right)\cap {\mathcal{R}}_{\mathrm{evader}}\left(t\right)\ne ⌀$.

The distinction becomes more significant under partial observability, where the adversary’s initial state ${\mathbf{x}}_a^0$ is also uncertain [69]. Wen and Gurfil [16] introduced the concept of RDs under initial state uncertainties, modeling ${\mathcal{X}}_0$ as an ellipsoid or polytope rather than a singleton. The resulting predictive set is the image of this initial uncertainty under dynamics, further bloated by unknown controls——a problem requiring set propagation techniques (zonotopes, polynomial approximations) rather than simple forward simulation.

The informational attribute thus introduces a crucial distinction often overlooked in classical RD literature: whether control inputs are known (deterministic planning) or unknown but bounded (predictive threat assessment). In orbital games, both perspectives are needed simultaneously: the ego spacecraft plans using its deterministic RD while estimating the adversary’s predictive RD to assess threats and identify opportunities. This dual perspective provides the theoretical foundation for the orbital game applications discussed in Section 4.

3. Solution Methods for Spacecraft Reachable Domain

Based on the unified modeling framework established in Section 2, this section reviews the core solution methods for computing spacecraft RD. Following the spatial-temporal classification introduced in Section 2, we distinguish between spatial scale and time scale RD methods. Within each category, we further categorize the methods according to their mathematical foundations. A comparative summary and critical analysis are provided at the end of this section to guide method selection for the orbital game applications analyzed in Section 4.

3.1. Methods for Spatial-Scale Reachable Domain

Spatial scale RD methods address the fundamental question: “What positions can the spacecraft eventually reach if given unlimited time?” These methods consider only fuel/energy constraints and are typically used in early mission design phases for capability assessment and target screening. Three families of methods have been developed, each with distinct mathematical foundations and computational characteristics.

3.1.1. Analytical and Semi-Analytical Methods

Analytical and semi-analytical methods exploit the integrability of two-body orbital dynamics to derive closed-form or semi-closed-form expressions for RD boundaries. They treat the problem as a geometric one——finding the envelope of a family of trajectories parameterized by control variables. Computationally efficient and geometrically insightful, these methods typically assume two-body dynamics and impulsive thrust, which are reasonable for many mission scenarios. Extensions to $J_2$ perturbed dynamics [70] and continuous-thrust approximations have also been developed, albeit with increased complexity.

Figure 3 illustrates the general framework shared by these methods. The common principle is to parameterize the RD boundary as an extremum problem over control variables such as maneuver point location and impulse direction.

The foundational work by Xue et al. [7] first classified the single-impulse RD into three scenarios: (i) fixed maneuver point with arbitrary impulse direction, (ii) arbitrary maneuver point with fixed impulse direction, and (iii) arbitrary maneuver point with arbitrary impulse direction. This classification established the basic framework for subsequent analytical developments. Three representative approaches have since been developed along this line.

The algebraic criterion method [8] transforms the reachability check into a sphere-hyperbola intersection problem. The set of all possible velocity vectors after an impulse is treated as a “maneuver sphere”, while the set of velocity vectors required to reach a target position is treated as a “terminal velocity hyperbola”. A target is reachable if and only if the hyperbola intersects or is tangent to the sphere. This geometric insight is converted into an algebraic accessibility function $\Delta \left(P\right)$, with $\Delta \left(P\right)\ge 0$ indicating reachability. The minimum required impulse magnitude is obtained by solving $\Delta \left(P\right)=0$ via Newton iteration. The method provides rapid feasibility assessment, but the numerical implementation is sensitive to initial guesses, and multiple solutions may exist without automatic criteria for selecting the physically relevant one.

The extremum envelope method [14] addresses the exact determination of RD boundaries. Rather than solving intersection conditions, it expresses the target position as the product of a direction vector and a distance, then establishes the reachable radius R as a function of control variables such as true anomaly f and impulse direction angles $\alpha$, $\beta$. For each fixed direction in space, the maximum reachable distance is found by solving the gradient-zero condition $\partial R/\partial \theta =0$ using Newton-Raphson iteration. Logarithmic coordinate transformation and neighborhood initial guess strategies are employed to ensure numerical stability. The envelope surface is constructed by collecting these extremal points across all directions. This method yields exact three-dimensional envelope surfaces, significantly outperforming the over-approximated results of Xue et al. [7], but still relies on numerical iteration for each boundary point.

The characteristic ellipse approximation method [15] completely avoids iterative computations by sacrificing a small amount of precision. The key insight is that the RD can be decomposed into two subdomains, with the inner and outer boundaries approximated by six characteristic ellipses, each given by closed-form expressions derived from orbital transfer geometry. The final RD is obtained by taking the convex hull of these six ellipses. This method achieves approximately 97% accuracy under typical conditions with extremely low computational cost, making it suitable for online real-time mission planning. The main limitation is that accuracy degrades when the impulse magnitude exceeds approximately 2000 m/s, as the underlying geometric assumptions break down.

Together, these three methods represent different trade-offs between geometric exactness and computational efficiency. The algebraic criterion is best suited for rapid feasibility checking, the extremum envelope for precise boundary determination requiring higher accuracy, and the characteristic ellipse approximation for real-time applications where speed is paramount.

3.1.2. Numerical Optimization Methods

Numerical optimization methods offer an alternative paradigm for RD computation. Unlike analytical methods that seek closed-form expressions through geometric insight, numerical optimization methods treat the RD problem as an optimal control problem, solving for trajectories that define the boundary. These methods are particularly valuable when the problem involves complex dynamics (e.g., $J_2$ perturbation), continuous thrust, or multiple constraints that resist analytical treatment. They offer greater flexibility at the cost of higher computational expense and dependence on numerical solvers.

Figure 4 shows the overall framework. The continuous optimal control problem is discretized or transformed and then solved using specialized algorithms. Early work by Wei and Li [10] employed numerical methods to investigate the absolute RD under continuous finite thrust, comparing the results with impulsive thrust cases and revealing significant morphological differences between the two thrust models. Two categories of numerical optimization methods are distinguished by their solution philosophy.

Indirect methods are based on Pontryagin’s maximum principle, which derives necessary conditions for optimality and transforms the problem into a two-point boundary value problem (BVP). Wang and Jiang [21] applied this approach to continuous low-thrust RDs within the framework of linearized relative motion equations (Tschauner-Hempel equations). Notably, their method yields analytical solutions for the RD envelope and its optimal control profile, demonstrating that indirect methods can, in certain linearized settings, produce closed-form expressions rather than purely numerical solutions. The Hamiltonian is formulated and optimality conditions are derived, yielding a BVP with boundary conditions specified at the initial and terminal times. The BVP is then solved numerically using shooting methods, typically requiring iterative adjustment of initial costate guesses.

In contrast, direct methods discretize the continuous-time optimal control problem and transform it into a nonlinear programming (NLP) problem, offering greater robustness to initial guesses. Using direct collocation methods, the time horizon is discretized (e.g., via Gauss pseudospectral or multiple shooting techniques), with the continuous dynamics approximated as algebraic constraints at the discretization nodes and the cost function approximated via quadrature. The resulting finite-dimensional NLP is then solved using standard solvers such as SNOPT or IPOPT. For reachable set computation, pseudospectral methods have been effectively employed for lunar landing reachability analysis [71] and cislunar low-thrust trajectory optimization [72]. The efficiency of these approaches can be further enhanced by exploiting sparsity in the resulting NLPs [73]. Direct methods are more robust to initial guesses than indirect methods because they do not require explicit derivation of optimality conditions or costate initialization, and they can readily handle complex constraints including state and control bounds, path constraints, and terminal conditions. The trade-off is that solution optimality depends on discretization resolution: coarse discretization may miss fine details of the RD boundary, while fine discretization leads to large NLP problems that are computationally expensive. Nevertheless, direct methods are widely used in engineering practice due to their flexibility and robustness.

3.1.3. Geometric and Sampling Methods

When analytical methods are inapplicable and numerical optimization methods are computationally prohibitive, sampling methods provide a pragmatic alternative. These methods treat the RD problem as a set estimation problem, approximating the reachable set through random or uniform sampling of the control space. They make few assumptions about the model but sacrifice exactness and efficiency.

Figure 5 presents the common workflow. The basic idea is to sample the control space, propagate the dynamics forward, and extract the RD boundary from the resulting point cloud.

The Monte Carlo sampling method is the most versatile approach. Control variables (impulse direction, magnitude, maneuver point, or thrust time history) are randomly or uniformly sampled according to their admissible ranges. The dynamics are propagated forward from the initial state for each sample, yielding a cloud of terminal state points. The RD boundary is then extracted from this point cloud using convex hull algorithms (for convex domains) or $\alpha$-shape algorithms (for non-convex domains). Du and Yang [74] used Monte Carlo simulation to verify the accuracy of their analytical model for single-impulse RDs. Sampling-based validation has also been employed in other studies to confirm analytical RD boundaries [9,14,75].

The main limitation of sampling methods is their slow convergence rate. For Monte Carlo sampling, the root mean square error of the estimate decreases as $O(1/\sqrt{N})$, where N is the number of samples. Achieving high precision therefore requires a large number of samples, which increases computational cost. Moreover, sampling methods provide no guarantees of covering the true RD—rare but important extremal trajectories may be missed by random sampling.

Nevertheless, their simplicity and minimal assumptions on the model make them attractive for validation purposes. Unlike analytical or numerical optimization methods, sampling methods do not require differentiability, convexity, or closed-form dynamics; they only require the ability to propagate the state forward for any given control input. This makes them particularly useful for verifying results obtained by other methods, especially when the problem involves complex dynamics or uncertainties.

3.1.4. Comparison and Selection Guide for Spatial-Domain Methods

The three families of spatial-domain methods presented above differ fundamentally in their mathematical principles, computational characteristics, and applicability domains. To provide an intuitive visual comparison, Figure 6 presents a radar chart evaluating these methods across five key dimensions: computational efficiency, solution accuracy, model flexibility, implementation complexity (lower is easier), and robustness under uncertainty.

The radar chart reveals the distinct trade-offs among the three method families. Analytical methods achieve the highest computational efficiency and provide valuable geometric insights, though their model flexibility is limited to two-body dynamics and impulsive thrust. Numerical optimization methods excel in solution accuracy and model flexibility, capable of handling complex dynamics and constraints, but incur the highest computational cost and implementation complexity. Sampling methods are distinguished by their dimension-independent convergence, which makes them well-suited for high-dimensional problems where other methods struggle with the curse of dimensionality. Their primary limitation is the slow convergence rate, requiring a large number of samples to achieve high accuracy. This makes them ideal for validation and exploratory analysis, though less suitable for high-precision or real-time applications.

Based on these characteristics, each method family is suited for specific application scenarios. Analytical methods are preferred for preliminary mission design and rapid capability assessment, where computational speed is prioritized over model fidelity. Numerical optimization methods are suitable for scenarios requiring high accuracy under complex dynamics, such as $J_2$ perturbation or continuous thrust. Sampling methods serve as validation tools and are valuable when the problem involves high-dimensional uncertainties or black-box dynamics that resist analytical treatment.

In practice, these methods can be used in a complementary manner. A potential workflow begins with analytical methods for initial screening to identify promising candidates. Numerical optimization can then be applied for precise boundary computation of the shortlisted candidates. Finally, sampling methods may validate the results under model uncertainties, providing additional confidence for mission planning decisions.

3.2. Methods for Time-Scale Reachable Domain

While spatial-scale methods focus on the maximum spatial envelope reachable under fuel or energy constraints, time-scale reachable domain (TSRD) methods address a more constrained and practically relevant problem: the set of positions that a spacecraft can reach at a specific time $t_f$ or within a specified time window $[t_1,t_2]$. Unlike free-time RDs, where the spacecraft can exploit arbitrarily long flight times to achieve energy-optimal transfers, TSRDs must satisfy both position and temporal constraints simultaneously. This coupling arises from Kepler’s equation, which relates transfer time to orbital geometry and admits no closed-form solution. Consequently, TSRDs are always proper subsets of their free-time counterparts:

$$\mathcal{R}\left(t_f\right)\subseteq {\mathcal{R}}_{\infty },\mathcal{R}\left([t_1,t_2]\right)=\bigcup _{\tau \in [t_1,t_2]}\mathcal{R}\left(\tau \right)\subseteq {\mathcal{R}}_{\infty }$$

(8)

Despite their importance in time-critical missions such as orbital interception, rendezvous at a prescribed epoch, and emergency return within a survivable window, the literature on TSRD remains relatively sparse. Most existing studies still rely on free-time approximations or linearized models that neglect the nonlinear coupling between time and position. To date, only a few dedicated TSRD solution methods have been developed. These methods can be classified into three categories according to their mathematical foundations and the trade-offs they strike between computational efficiency and model fidelity.

3.2.1. Linearized Ellipsoidal Approximation Method

For scenarios involving small impulse magnitudes and short propagation times, Ma and Zhang [31] proposed an ellipsoidal approximation method based on the linearized state transition matrix (STM). The key idea is to linearize the two-body dynamics around a nominal trajectory, which allows the impulsive RD to be approximated as a linear transformation of the initial maneuver sphere.

Under this linearization, the deviation vector $\Delta {\mathbf{r}}_f$ at terminal time $t_f$ is related to the impulse vector $\Delta \mathbf{v}$ by:

$$\Delta {\mathbf{r}}_f={\Phi }_{12}\Delta \mathbf{v}$$

(9)

where ${\Phi }_{12}\in {\mathbb{R}}^{3\times 3}$ is the position-velocity sub-block of the two-body STM. The set of admissible impulses forms a sphere: $\parallel \Delta \mathbf{v}\parallel \le \Delta v_{max}$. Under linear mapping, this sphere becomes an ellipsoid in position space:

$$\Delta {\mathbf{r}}_f^{\mathsf{T}}\Lambda \Delta {\mathbf{r}}_f=1,\Lambda ={\left({\Phi }_{12}^{-1}\right)}^{\mathsf{T}}\left({\displaystyle \frac{1}{\Delta v_{max}^2}}I\right){\Phi }_{12}^{-1}$$

(10)

This ellipsoid serves as an approximate description of the TSRD envelope at time $t_f$. The method is computationally efficient, requiring only matrix multiplications and a simple root-finding step, making it suitable for real-time applications such as onboard pursuit impulse calculation. However, its accuracy is fundamentally limited by the linearization assumption. The linearized approximation is valid only when the propagation time is sufficiently short relative to the orbital period.

3.2.2. Exact Envelope Determination Method

To overcome the limitations of linearized approximations, Zhang et al. [29] developed an exact method for TSRD computation that does not rely on linearization. This method treats the TSRD envelope as the solution to a constrained extremum problem in the space of maneuver parameters.

The approach begins by parameterizing the transfer orbit orientation using two angles: $\beta$ (the rotation angle between the initial orbital plane and the transfer plane) and $\theta$ (the variation of true anomaly along the transfer orbit). For a given $\beta$, the admissible range of $\theta$ at time $t_f$ is determined by solving Kepler’s equation via an unconstrained nonlinear programming formulation:

$$min\left\{\Delta f_{min}\left(\alpha \right)=\Delta f,\Delta f_{max}\left(\alpha \right)=-\Delta f\right\}$$

(11)

where $\alpha$ is the angle describing the direction of the impulse within the transfer plane. The extremal values of the terminal radius $r_f$ for each orientation $(\beta ,\theta )$ are then obtained by solving the time constraint equation $\Delta t\left(\alpha \right)=t_f$ using numerical root-finding (e.g., Newton-Raphson iteration). The complete TSRD envelope is constructed by collecting these extremal points across all feasible orientations.

The method distinguishes between two subcategories of TSRD: the RD at a given time $\mathcal{R}\left(t_f\right)$ and the RD over a duration $\mathcal{R}\left([t_1,t_2]\right)$. For the latter, the equality constraint $\Delta t=t_f$ is replaced by inequality constraints $t_1\le \Delta t\le t_2$, and the envelope is obtained by solving a constrained nonlinear programming problem.

The exact method achieves high accuracy, with the computed envelope closely matching Monte Carlo simulation results. However, this accuracy comes at a significant computational cost: solving a single TSRD envelope requires 3–4 seconds on a modern processor, and the algorithm involves nested iterations over $\beta$, $\theta$, and $\alpha$. This computational burden limits its applicability to offline analysis or scenarios where real-time performance is not required.

3.2.3. Fast Analytical Approximation via Series Expansion

Recognizing the need for a method that balances accuracy and efficiency, Yang et al. [30] proposed a fast analytical approximation for TSRD computation. This method leverages the universal variable formulation of Lambert’s problem and a second-order series expansion of Kepler’s equation, thereby avoiding numerical iterations altogether.

The fundamental innovation lies in the analytical treatment of the eccentric anomaly E. Using the series expansion method valid for orbits with eccentricity $e<0.5$, the eccentric anomaly is expressed as:

$$E=M+\left(e-{\displaystyle \frac{e^3}{8}}\right)sinM+{\displaystyle \frac{e^2}{2}}sin2M+{\displaystyle \frac{3e^3}{8}}sin3M+O\left(e^4\right)$$

(12)

where $M=M_0+n\Delta t$ is the mean anomaly. Substituting this expansion into the universal variable formulation yields closed-form expressions for the Lagrange coefficients f and g, which then allow the terminal position to be computed as:

$$\mathbf{r}\left(t_f\right)=f(t_f,\Delta \mathbf{v}){\mathbf{r}}_0+g(t_f,\Delta \mathbf{v})({\mathbf{v}}_0+\Delta \mathbf{v})$$

(13)

For a given flight time $t_f$ and impulse magnitude $\parallel \Delta \mathbf{v}\parallel$, the terminal position $\mathbf{r}\left(t_f\right)$ becomes an explicit function of the impulse direction angles $(\alpha ,\beta )$. The TSRD envelope is then obtained by uniformly sampling the unit sphere of impulse directions, a process that takes only milliseconds to complete.

The method is particularly well-suited for near-circular orbits and short-to-medium propagation times. Under these conditions, the position error introduced by the second-order truncation remains relatively modest, and is generally not the dominant source of trajectory deviation compared to the effects of orbital perturbations and measurement uncertainties. Moreover, any residual truncation error can be compensated by a rolling feedback mechanism when the method is embedded in a closed-loop guidance framework. For orbits with larger eccentricities, the series expansion can be extended to higher orders to maintain accuracy, albeit with a moderate increase in computational cost. This tunable trade-off between accuracy and speed makes the fast analytical approximation method a practical choice for real-time onboard applications.

3.2.4. Comparative Summary of TSRD Methods

To illustrate the applicability of the linearized ellipsoidal approximation, Figure 7 compares the TSRD envelopes obtained by the STM-based linearized method (Method 1) and the fast analytical approximation (Method 2) at short and long propagation times, with Monte Carlo simulation results serving as the ground truth for reference.

As shown in Figure 7, at short propagation times, both methods yield comparable results. However, when the propagation time increases, the STM-based method exhibits obvious deviation, whereas the fast analytical approximation retains its accuracy. This confirms that the STM-based linearized approximation is suitable only for short-time applications.

Regarding computational efficiency, the three methods exhibit distinct characteristics. The linearized ellipsoidal approximation [31] is the most efficient, as it relies solely on matrix operations and a simple root-finding step. The fast analytical approximation [30] also achieves low computational cost by avoiding numerical iterations through series expansion, making it suitable for real-time onboard applications. The exact envelope determination method [29], while offering the highest accuracy, requires nested numerical iterations and solving Kepler’s equation, resulting in significantly higher computational expense.

Given the current scarcity of TSRD literature, these three methods represent the state of the art. In practical applications, a hybrid strategy can be adopted: the STM-based method is recommended for short-time scenarios where efficiency is prioritized, while the fast analytical approximation or the exact envelope determination should be employed for long-time scenarios where accuracy becomes critical.

4. Applications of Reachable Domain in Orbital Games

4.1. Fundamentals of Orbital Games and the Role of Reachable Domain

4.1.1. Typical Scenarios and Classification of Orbital Games

Orbital games refer to scenarios in which two or more spacecraft engage in competitive or cooperative interactions through orbital maneuvers, typically involving pursuit, evasion, encirclement, or protection. These problems have attracted increasing attention due to the growing demand for on-orbit servicing, active debris removal, and space situational awareness.

Depending on the number of participating spacecraft and their interaction patterns, orbital games can be classified into three typical categories.

The first category is pursuit-evasion games, which represent the most fundamental orbital game paradigm. In this setting, a pursuer spacecraft aims to intercept or rendezvous with an evader, while the evader attempts to escape or avoid capture. This problem is closely related to interception feasibility, fuel consumption, and time-critical mission constraints [31,32].

The second category is cooperative encirclement, where multiple chasing spacecraft coordinate their maneuvers to surround or capture a non-cooperative target. This scenario is particularly relevant to multi-agent space missions, such as cooperative inspection, coordinated interception, and formation-based capture [38,67].

The third category is threat avoidance and defense, in which an ego spacecraft must evade potential threats posed by an adversarial or non-cooperative object. The objective is to maintain a safe separation distance or to avoid entering regions where the adversary’s capture capability is effective [43,68].

These three categories share a common mathematical core: in each case, the spacecraft’s maneuvering capability is bounded by fuel or thrust constraints, and the interaction outcome is determined by whether one spacecraft’s reachable set intersects with another’s trajectory or reachable set. This inherent connection provides a unifying perspective for analyzing orbital games through the lens of RD theory.

The three scenarios illustrated in Figure 8 demand different yet interconnected uses of RD theory. In pursuit-evasion games, the core question is whether the pursuer’s RD intersects the evader’s reachable set within the available time and fuel budget. In cooperative encirclement, the focus shifts to how multiple pursuers’ RDs can collectively cover the target’s escape routes. In threat avoidance, the emphasis is on keeping the ego spacecraft’s trajectory disjoint from the adversary’s predictive RD. The following sections will examine, in turn, how RD methods have been applied to each of these scenarios.

4.1.2. Dual Role of Reachable Domain in Orbital Games

As established in Section 2.4, the RD assumes two fundamentally different roles depending on whose control inputs are being considered: the deterministic RD for known maneuvers and the predictive RD for unknown but bounded adversarial maneuvers. In orbital game scenarios, this distinction translates directly into a dual-use framework.

On one hand, each participating spacecraft can compute its own deterministic RD to support proactive decision-making. For a pursuer, this domain answers the question: "Given my fuel budget and thrust capability, what target states can I reach?" For an evader, it informs: "What escape trajectories are available to me?" This perspective corresponds to planning under complete information about one’s own control [22,76].

On the other hand, when confronting a non-cooperative adversary with unknown maneuver timing or direction, the predictive reachable domain serves as a threat envelope. The pursuer can estimate the set of all possible states the evader could reach using only the bound on its maneuvering capability. This perspective corresponds to threat assessment under partial information [69,77].

4.2. One-on-One Pursuit-Evasion Games

The fundamental problem in one-on-one pursuit-evasion games is determining whether a pursuer can intercept an evader given their respective maneuvering capabilities. The RD framework provides an elegant answer: interception is feasible if and only if the pursuer’s RD intersects the evader’s reachable set within the mission constraints. This set-theoretic formulation transforms a dynamic game problem into a geometric condition, enabling systematic analysis of interception feasibility, fuel requirements, and time constraints. More specifically, a sufficient condition for the pursuer to guarantee interception is that the evader’s entire RD is contained within that of the pursuer. Conversely, if the pursuer’s RD is entirely disjoint from the evader’s reachable set, successful interception is impossible. These conditions can be expressed formally as:

$$\mathrm{Guaranteed}\mathrm{interception}\text{:}{\mathcal{R}}_e\left(t\right)\subseteq {\mathcal{R}}_p\left(t\right),\forall t\in [t_0,t_f],$$

$$\mathrm{Impossible}\mathrm{interception}\text{:}{\mathcal{R}}_p\left(t\right)\cap {\mathcal{R}}_e\left(t\right)=\varnothing ,\forall t\in [t_0,t_f],$$

where ${\mathcal{R}}_p\left(t\right)$ and ${\mathcal{R}}_e\left(t\right)$ denote the RDs of the pursuer and evader at time t, respectively. Between these two extremes lies the regime of probabilistic or conditional interception, where the outcome depends on the specific maneuvers chosen by the evader.

For free-time interception where no time constraint is imposed, the condition reduces to whether the union of the pursuer’s RDs over all future times overlaps with that of the evader. Ma and Zhang [31] quantified this condition through RD coverage analysis: the minimum Delta-V required for interception is determined by the geometric overlap between two ellipsoidal RDs under the Clohessy-Wiltshire framework. Wang et al. [34] extended this line of research by developing a feasibility analysis method for orbital pursuit-evasion missions based on single-impulse RDs, providing practical criteria for mission planning. For continuous thrust scenarios, Gong et al. [32] established that interception is guaranteed if the pursuer’s RD eventually encloses the evader’s entire trajectory, which is a dynamic extension of the set containment condition stated above. Bowerfind and Taheri [23] developed rapid approximation methods for low-thrust reachable sets using indirect multi-stage algorithms, achieving near-real-time computation suitable for onboard interception feasibility assessment.

When interception must occur at a specific time or within a time window, the fixed-time RD becomes essential, and the containment condition must be evaluated at the prescribed terminal time rather than over all future times. Zhang et al. [29] systematically distinguished two subproblems: interception at a prescribed time $t_f$ and interception within a window $[t_1,t_2]$. For the former, the condition reduces to a membership test of the evader’s predicted position in ${\mathcal{R}}_p\left(t_f\right)$. For the latter, one must evaluate whether ${\mathcal{R}}_e\left(\tau \right)\subseteq {\mathcal{R}}_p\left(\tau \right)$ for some $\tau \in [t_1,t_2]$, which is computationally more demanding as it requires analyzing the union of RDs over a continuous interval. The authors developed an exact envelope determination method that solves Kepler’s equation numerically, providing high accuracy for offline analysis. For scenarios requiring real-time decision-making, Zhang et al. [36] proposed an analytical impulse thrust strategy for spacecraft pursuit-evasion games that accounts for mass variation and perception frequency. Their method leverages the RD to compute optimal impulse timing and direction without iterative numerical optimization. From the learning perspective, Qiao et al. [41] introduced a time-limited RD-guided learning method for spacecraft game decision-making. This approach integrates RD computations with deep reinforcement learning: the RD provides a safety certificate that guides exploration during training, while the neural network approximates the optimal policy for real-time execution. The time-limited RD is particularly well-suited for this purpose because it respects mission time constraints, which are often critical in orbital pursuit-evasion scenarios.

Beyond the basic interception problem, the RD framework also supports the analysis of equilibrium strategies in pursuit-evasion games. Li et al. [78] developed a Nash equilibrium solution method based on the relative motion reachable set. Their key contribution is transforming the pursuit-evasion game into a set-based optimization problem where equilibrium corresponds to mutual RD conditions. Xu et al. [79] extended this approach to continuous thrust scenarios, proposing numerical methods for Nash equilibrium using RD analysis. The non-convexity of continuous thrust RDs requires global optimization techniques, making the problem more computationally demanding but also more realistic for high-fidelity mission analysis. A critical engineering decision in any such analysis is the choice between impulsive and continuous thrust models. Impulsive models, such as those analyzed by Ma and Zhang [31] and Wang et al. [34], assume instantaneous velocity changes and produce convex RDs, enabling efficient interception condition checks including the containment test ${\mathcal{R}}_e\subseteq {\mathcal{R}}_p$. Continuous thrust models, as studied by Gong et al. [32] and Bowerfind and Taheri [23], yield non-convex RDs but capture the realistic evolution of thrust over time. For practical mission design, the recommended strategy is to use impulsive models for preliminary feasibility assessment and real-time guidance, while reserving continuous thrust models for high-fidelity offline analysis where accuracy is paramount.

4.3. Multi-Agent Cooperative Games

In contrast to one-on-one pursuit-evasion scenarios, multi-agent cooperative games involve multiple spacecraft working together to achieve a common objective, such as surrounding a non-cooperative target or coordinating interception from multiple directions. The RD framework extends naturally to these settings by considering the collective coverage of multiple pursuers’ reachable sets. A fundamental condition for successful cooperative encirclement is that the target’s entire RD is covered by the union of all chasers’ RDs:

$${\mathcal{R}}_t\left(t\right)\subseteq {\displaystyle \bigcup _{i=1}^n}{\mathcal{R}}_i\left(t\right),\forall t\in [t_0,t_f],$$

where ${\mathcal{R}}_t\left(t\right)$ is the target’s RD and ${\mathcal{R}}_i\left(t\right)$ is the RD of the i-th chaser. When this condition holds, the target has no escape opportunity regardless of its maneuver choices, representing a cooperative analog of the guaranteed interception condition in one-on-one games. Conversely, if a gap exists in the coverage, the target may exploit that gap to evade capture.

The design of cooperative encirclement configurations has been investigated by several researchers. Li et al. [37] studied cooperative encirclement of non-cooperative targets using spacecraft swarms, proposing a semi-analytical solution for the cooperative RD. Their analysis revealed that three chasers forming an equilateral triangular configuration provide optimal coverage for planar encirclement scenarios, as this arrangement minimizes the uncovered region while using the minimum number of chasers. The coverage condition can be verified by checking whether the target’s reachable ellipsoid is completely contained within the Minkowski sum of the chasers’ reachable ellipsoids. For scenarios requiring higher reliability, additional chasers can be deployed to create overlapping coverage regions, providing redundancy against individual chaser failures or unexpected target maneuvers.

When the target is not stationary but actively maneuvering, the encirclement problem becomes more challenging because the coverage condition must hold over an extended time horizon. Li and Zhang [38] extended the analysis to relative RDs in target-centered reference frames, computing for a given target orbital state the set of relative positions from which a chaser can intercept any potential escape trajectory. This "capture region" serves as a guidance map for positioning chasers to maximize encirclement effectiveness. The authors demonstrated that by maintaining each chaser’s position within its respective capture region, the coalition can guarantee that the target remains encircled despite its maneuvering attempts. Lin and Zhang [25] further extended this analysis to elliptical reference orbits using continuous-thrust reachable sets for linear relative motion described by the Tschauner-Hempel equations, providing more accurate encirclement guidance for missions in highly elliptical orbits where the circular-orbit approximation is inadequate.

Beyond encirclement, the RD framework also supports cooperative interception where multiple interceptors are tasked with capturing one or more targets. Qing and Zhang [39] studied multi-to-multi spacecraft impulsive cooperative interception based on RD coverage. Their task allocation formulation treats each target as covered if at least one interceptor’s RD intersects the target’s reachable set within the mission time window. The overall mission success probability can be estimated from individual interception probabilities:

$$P_{\mathrm{success}}=1-{\displaystyle \prod _{i=1}^n}\left(1-p_i\right),$$

where $p_i$ is the interception probability of the i-th interceptor, estimated from the volumetric overlap between the interceptor’s RD and the target’s reachable set. This formulation assumes independence of interception attempts, which is a reasonable approximation when interceptors operate from distinct initial positions and the target’s maneuver strategy is unknown. For scenarios with a single highly maneuverable target, Wang and Zhang [40] developed an impulsive maneuver strategy that uses RD predictions to guide interceptor positioning without requiring explicit target maneuver models, achieving robust performance under sparse reward conditions through a concept they term "RD compression."

For multi-pursuer/one-evader scenarios, Li et al. [80] proposed a geometric method based on RD theory. Their approach transforms the Nash equilibrium condition into a RD equivalence representation, where the equilibrium terminal time is determined by checking the geometric relationship between the pursuers’ reachable sets and the evader’s reachable set. Using grid-point search to compute RD boundaries and solid angle calculation to quantify coverage, the method reduces computational time from over two hours required by traditional differential game methods to approximately twelve minutes, demonstrating the practical advantage of reachable set-based approaches.

A critical consideration in multi-agent cooperative games is the coordination mechanism used to achieve collective coverage. Centralized coordination, where a single agent computes the optimal coverage configuration and assigns positions to all chasers, is conceptually straightforward and can achieve global optimality [37,39]. However, it suffers from single-point-of-failure vulnerability and requires reliable communication links. Decentralized coordination, where each spacecraft computes its own RD and negotiates coverage through local information exchange, offers greater robustness and scalability at the cost of optimality guarantees [40]. Hybrid approaches that combine offline centralized planning with online decentralized adjustment have also been explored, offering a pragmatic trade-off between optimality and robustness for practical mission design.

4.4. Threat Avoidance and Defense Games

While Section 4.2 and Section 4.3 focus on the pursuer’s perspective of intercepting or encircling a target, this section addresses the dual problem from the evader’s or defender’s perspective: how to avoid capture, evade threats, or protect high-value assets from adversarial approach. The earliest work recognizing that reachable domains can serve as a tool for safety analysis dates back to the early 2010s [42], establishing the foundation for subsequent threat avoidance studies. Threat avoidance and defense games are particularly relevant for space asset protection, on-orbit servicing safety, autonomous collision avoidance, and spacecraft operation in debris-populated environments. The fundamental condition for safe operation is that the ego spacecraft’s trajectory remains disjoint from the adversary’s RD, augmented by a safety buffer to account for navigation errors and control uncertainties:

$${\mathbf{x}}_{\mathrm{ego}}\left(t\right)\notin {\mathcal{R}}_{\mathrm{adv}}\left(t\right)\oplus {\mathcal{B}}_{\mathrm{safety}},\forall t\in [t_0,t_f],$$

where ${\mathcal{R}}_{\mathrm{adv}}\left(t\right)$ is the adversary’s predictive RD and ${\mathcal{B}}_{\mathrm{safety}}$ is a safety buffer zone (e.g., a sphere of radius $r_{\mathrm{safe}}$). Conversely, an adversary is considered successfully deterred if its predictive RD has no intersection with the asset’s RD: ${\mathcal{R}}_{\mathrm{adv}}\left(t\right)\cap {\mathcal{R}}_{\mathrm{asset}}\left(t\right)=\varnothing$. Between these two extremes lies a gray zone where the outcome depends on maneuver timing, information asymmetry, or relative fuel budgets.

The core idea of using RDs for threat avoidance is to treat the adversary’s maneuvering capability as a threat envelope that the ego spacecraft must avoid. Zhang et al. [35] proposed a threat avoidance strategy for spacecraft maneuvering approach based on orbital RD analysis, demonstrating that a safe approach trajectory can be incrementally constructed by ensuring the ego spacecraft remains outside the target’s predictive RD at all times. This predictive reachable domain is computed based on the target’s maximum maneuvering capability, which is typically the only available information in non-cooperative settings. The method produces safer trajectories than traditional potential field methods, particularly when the target’s maneuver timing is uncertain. Building on this work, Zhang et al. [43] extended the analysis to orbital threat assessment, providing methods to compute the minimum safe distance and identify safe time windows during which the threat level is acceptably low. These safe windows are particularly useful for mission planning in cluttered or adversarial environments, as they allow the spacecraft to schedule critical operations during periods of reduced risk.

The first extension considers passive threats, where the adversary is not actively maneuvering but poses a collision risk due to its uncontrolled motion. Wen and Gurfil [44] applied the RD method to model the early-to-medium-term evolution of debris clouds following satellite breakups. By treating each debris fragment’s RD as an expanding set, they demonstrated that the overall debris cloud’s spatial distribution could be efficiently approximated without propagating individual trajectories. This approach enables spacecraft operators to assess collision risk from debris clouds and plan avoidance maneuvers accordingly. Xu et al. [45] proposed a zonotope-based method for computing reachable sets of satellite relative motion, achieving a balance between computational efficiency and geometric accuracy. Zonotopes preserve convexity while capturing more geometric detail than ellipsoidal approximations, enabling real-time collision checking for proximity operations where safety margins are measured in meters rather than kilometers. For collision avoidance against debris objects, Li and Zhang [46,47] proposed both continuous-thrust and impulsive optimal avoidance maneuvers based on RD analysis. Their key result is that the minimum required Delta-V for collision avoidance can be computed by solving for the shortest distance between the ego spacecraft’s RD and the threat object’s predicted trajectory. This formulation unifies collision avoidance with interception analysis: interception seeks to reduce the distance between RDs to zero, while avoidance seeks to maximize that distance.

The second extension addresses active threats, where the adversary maneuvers intelligently to pursue or intercept the ego spacecraft. From the evader’s perspective, this becomes a pursuit-evasion game where the goal is to maximize the distance from the pursuer’s RD. Zhang et al. [48] developed an escape-zone-based optimal evasion guidance strategy against multiple orbital pursuers. Their key contribution is the concept of the "escape zone"—the set of states from which the evader can successfully avoid capture regardless of the pursuers’ strategies. By computing the boundary of this escape zone using RD analysis, the evader can identify safe regions in the state space and plan trajectories that maximize the minimum distance to the pursuers’ RDs.

The third extension considers cooperative defense, where a single spacecraft’s avoidance capability is insufficient and multiple guardians must coordinate to protect a high-value asset. Zhang et al. [49] proposed a cooperative guard strategy based on RD coverage. The guardians maintain a defense perimeter by positioning themselves such that the union of their RDs completely covers the asset’s vulnerable directions. The required number of guardians scales linearly with the angular span of the asset’s vulnerable directions, and the authors derived conditions for minimal guardian configuration. When the adversary’s predictive RD can be estimated, the protection condition becomes ${\mathcal{R}}_{\mathrm{adv}}\left(t\right)\cap {\mathcal{R}}_{\mathrm{asset}}\left(t\right)=\varnothing$, which can be maintained by the guardians through continuous adjustment of their positions to close emerging coverage gaps. This cooperative defense paradigm is particularly relevant for protecting space stations, critical communication satellites, or on-orbit servicing facilities.

The fourth extension incorporates real-world dynamical effects that alter the shape and orientation of threat envelopes. Gao et al. [81] addressed threat region assessment under $J_2$ perturbation, comparing the RD distribution under the two-body model and the $J_2$-perturbed model. Their analysis revealed that Earth’s oblateness significantly affects the spatial distribution of the RD, altering both its size and orientation. Based on this finding, they further investigated the spatial rendezvous characteristics between a non-cooperative target’s RD and the orbits of potential threats, enabling quantitative threat region assessment. To comprehensively evaluate medium-to-long-term collision risks, they constructed an assessment index system that includes Effective Coverage Rate, Distribution Concentration, and Distribution Skewness. This framework provides quantitative support for space situational awareness and orbit protection decision-making in low Earth orbit where $J_2$ effects cannot be neglected.

The methods reviewed in this section cover passive threats, active evasion, cooperative defense, and perturbation effects. Together, they provide a set of practical approaches for designing safe spacecraft operations in adversarial or uncertain environments. Combined with the pursuit and encirclement methods from Section 4.2 and 4.3, these approaches complete a unified RD framework for orbital game analysis from both offensive and defensive perspectives.

5. Challenges and Future Directions

5.1. Limitations of Existing Approaches

Despite significant progress in RD modeling, computation, and applications, several limitations persist that constrain the applicability of existing methods to real-world orbital game scenarios. These limitations span three interconnected aspects: modeling fidelity, computational efficiency, and scalability under uncertainty.

The first limitation concerns the modeling fidelity of existing RD methods. The vast majority of studies assume two-body dynamics or linearized relative motion, with only preliminary extensions to $J_2$ perturbation [70] and the circular restricted three-body problem [56]. Most methods further assume impulsive thrust [8,14] or simplified continuous-thrust models [32], neglecting the practical constraints of finite specific impulse, thrust magnitude limits, and attitude control coupling. These modeling simplifications, while enabling analytical tractability, limit the applicability of RD methods to real-world orbital scenarios where perturbations, nonlinearities, and operational constraints cannot be ignored.

The second limitation is the persistent trade-off between computational efficiency and solution accuracy. High-fidelity numerical optimization methods can handle complex dynamics and constraints but require substantial computation, rendering them unsuitable for onboard real-time applications. Fast analytical approximations, such as the characteristic ellipse method [15] and linearized ellipsoidal approximations, achieve millisecond-level computation but their accuracy degrades under long propagation times, large maneuvers, or highly nonlinear dynamics. Sampling methods, while flexible and dimension-independent, suffer from slow convergence and provide no guarantees of covering extremal trajectories. This trade-off remains particularly acute for time-critical orbital games, where both accuracy and speed are essential.

The third limitation concerns the scalability of RD methods to multi-agent settings and their ability to handle comprehensive uncertainties. Existing studies on cooperative encirclement and interception are largely restricted to small teams of three to five spacecraft [37,39], with limited investigation of large-scale swarms or distributed coordination. Furthermore, current approaches handle initial state uncertainty [16], unknown adversarial maneuvers [69], and dynamical perturbations [70] separately or under restrictive assumptions. A unified framework that simultaneously accounts for these multiple sources of uncertainty is lacking, limiting the practical deployment of RD methods for robust orbital game analysis under realistic operating conditions.

5.2. Future Research Directions

To address the limitations outlined above, several research directions warrant further investigation. These directions span complex dynamical environments, real-time computation, multi-agent scalability with initial configuration optimization, uncertainty handling, and experimental validation.

A pressing need is the extension of RD methods to complex dynamical environments. Building upon recent efforts [56], future work should develop computational frameworks for RDs in the circular restricted three-body problem. The invariant manifolds associated with Lagrangian points offer a natural structure for characterizing reachable sets in cislunar space, with potential applications to low-energy transfer design, orbit protection, and threat assessment near the Earth-Moon L1 and L2 points. Beyond CR3BP, other perturbation regimes such as high-fidelity ephemeris models and solar radiation pressure also deserve attention.

Equally important is the development of real-time reachable domain computation methods. Machine learning techniques offer promising pathways in this direction. Physics-informed neural networks can embed known dynamical constraints into the learning process, reducing data requirements and improving generalization. Reduced-order surrogate models, trained on offline optimal control solutions, can approximate RD boundaries in milliseconds, bridging the gap between numerical accuracy and onboard efficiency. Hybrid approaches that combine analytical approximations [15] with learning-based corrections may offer the best of both worlds, achieving accuracy comparable to numerical methods at a fraction of the computational cost.

Beyond computational challenges, the scalability of reachable domain methods to multi-agent settings remains a critical issue. Distributed RD computation and consensus algorithms that scale to large spacecraft swarms require further development, as existing centralized coordination mechanisms become computationally prohibitive when the number of agents grows. A specific problem within this broader scalability challenge is the initial configuration optimization for pursuit spacecraft based on RD coverage. Existing studies [37,39] have demonstrated that three chasers in triangular formation can achieve effective coverage, but systematic methods for optimizing the number, initial positions, and relative geometry of pursuers given an arbitrary target RD have not been established. This inverse problem has direct implications for cooperative encirclement, constellation deployment, defensive formation design, and multi-target interception. Future research could formulate this as a constrained optimization problem, with objective functions such as maximizing covered volume, minimizing the number of pursuers, or achieving uniform coverage subject to fuel and communication constraints.

Another important direction concerns the handling of uncertainties. Current approaches treat initial state uncertainty [16], unknown adversarial maneuvers [69], and dynamical perturbations [70] separately or under restrictive assumptions. A unified framework that integrates these multiple sources of uncertainty within a single probabilistic or set-theoretic formulation is needed. Robust RD methods that provide guaranteed performance bounds under worst-case uncertainties would be particularly valuable for safety-critical applications such as collision avoidance and asset protection. Stochastic RDs, which characterize the probability distribution of terminal states rather than just the set of possible states, offer a natural framework for risk-aware decision-making. Extensions to non-parametric uncertainty models and adversarial uncertainty representations, where the evader’s strategy is intentionally deceptive, represent challenging but important directions.

Finally, the intersection of RD theory with learning-based decision-making and experimental validation offers a promising pathway toward practical deployment. RDs can provide verifiable safety guarantees and reward shaping for deep reinforcement learning, guiding exploration toward feasible and safe regions of the state space while penalizing unsafe actions. Preliminary work [41] has demonstrated the feasibility of RD-guided learning for spacecraft game decision-making, but broader integration remains to be explored. Conversely, learning could accelerate RD computation by amortizing expensive trajectory optimization across similar problem instances, enabling meta-learning of RD boundaries that adapt to new scenarios with minimal computation. On the experimental front, hardware-in-the-loop demonstrations and on-orbit experiments are needed to validate RD computations under realistic conditions, including sensor noise, communication delays, actuator imperfections, and model errors. Such validation would accelerate the transition of reachable domain methods from academic research to operational space systems.

6. Conclusions

This paper has presented a comprehensive review of the spacecraft RD, encompassing its unified modeling framework, solution methods, and applications in orbital games. The review follows a coherent progression from fundamental theory to computational methods and finally to practical applications, providing a systematic reference for researchers in astrodynamics and space situational awareness.

A unified mathematical framework for the spacecraft reachable domain is established through three complementary classification dimensions. The spatial dimension distinguishes absolute reachable domains, which characterize overall maneuvering capability in inertial space, from relative RDs, which focus on position sets with respect to a target. The temporal dimension differentiates free-time RDs, representing the maximum spatial envelope under fuel constraints, from fixed-time RDs, which satisfy both position and time constraints for time-critical missions. The informational dimension distinguishes deterministic RDs, where control inputs are fully known, from predictive RDs, where only bounds on adversarial maneuvers are available. This three-dimensional taxonomy provides a systematic framework for selecting appropriate modeling approaches across diverse mission scenarios.

The solution methods for computing spacecraft RDs are organized according to the spatial-temporal classification. For spatial-domain RDs, three families of methods are examined. Analytical and semi-analytical methods, including the algebraic criterion method, extremum envelope method, and characteristic ellipse approximation, offer closed-form or semi-closed-form solutions with high computational efficiency, making them suitable for preliminary mission design and real-time applications. Numerical optimization methods, both indirect and direct, provide flexibility for complex dynamics and constraints at the cost of higher computation, making them suitable for high-fidelity offline analysis. Geometric and sampling methods, including Monte Carlo approaches, make minimal assumptions about the model and serve as valuable validation tools. For time-scale RDs, three dedicated methods are analyzed: the linearized ellipsoidal approximation based on the state transition matrix, which is efficient but limited to short propagation times; the exact envelope determination method solving Kepler’s equation numerically, which achieves high accuracy at significant computational cost; and the fast analytical approximation using series expansion, which balances accuracy and efficiency for near-circular orbits and short-to-medium propagation times.

The applications of RD theory to orbital games are examined through three representative scenarios. For one-on-one pursuit-evasion games, the fundamental interception condition is expressed as set intersection, with guaranteed interception requiring containment of the evader’s RD within that of the pursuer, and guaranteed safety requiring disjointness. Time-constrained interception is addressed through fixed-time RD analysis, with distinctions made between prescribed-time and time-window interception. For multi-agent cooperative games, cooperative encirclement requires the target’s RD to be covered by the union of all chasers’ RDs, with triangular configurations providing optimal planar coverage. Cooperative interception extends this framework to task allocation and success probability estimation. For threat avoidance and defense games, four extensions to the basic avoidance framework are identified: passive threats (debris clouds and collision avoidance), active threats (escape zone guidance and optimal evasion), cooperative defense (multi-guardian asset protection), and perturbation effects ($J_2$-influenced threat assessment). The duality between pursuit and avoidance is formalized through RD containment: pursuers seek ${\mathcal{R}}_{\mathrm{evader}}\subseteq {\mathcal{R}}_{\mathrm{pursuer}}$ for interception, while evaders seek ${\mathcal{R}}_{\mathrm{evader}}\cap {\mathcal{R}}_{\mathrm{pursuer}}=\varnothing$ for safety.

This review identifies several limitations of existing approaches, including modeling fidelity (over-reliance on two-body dynamics and simplified thrust models), computational trade-offs (accuracy versus speed), and scalability issues (multi-agent complexity and fragmented uncertainty handling). To address these limitations, future research directions are proposed, including extensions to complex dynamical environments such as the circular restricted three-body problem, real-time computation through machine learning, initial configuration optimization for cooperative pursuit, unified frameworks for multi-source uncertainty, and experimental validation through hardware-in-the-loop and on-orbit demonstrations.

Together, the unified modeling taxonomy, solution methods, and application scenarios reviewed in this paper establish a solid foundation for RD research. Addressing the challenges outlined above will advance the state of the art in space situational awareness, autonomous rendezvous, and orbital security.