Escape-Zone-Based Optimal Evasion Guidance Against Multiple Orbital Pursuers

The orbital evasion problem is getting increasing attention because of the increase of space maneuvering objects. In this article, an escape-zone-based optimal orbital evasion guidance law for an evading spacecraft on near circular reference orbit is proposed against multiple pursuing spacecraft with impulsive thrust. The relative reachable domain is introduced first and approximated as an ellipsoid propagating along the nominal trajectory under the short-term assumption. The escape zone for the impulsive evasion problem is presented herein as a geometric description of the set of terminal positions for all the impulsive evasion trajectories that are not threatened by the maneuvers of pursuers at the maneuver moment. A general method is developed next to calculate the defined escape zone through finding the intersection of two relative reachable domain approximate ellipsoids at arbitrary intersection moment. Then, the two-sided optimal strategies for the orbital evasion problem are analyzed according to whether the escape zone exists, based on which the escape value is defined and used as the basis of the proposed orbital evasion guidance scheme. Finally, numerical examples demonstrate the usefulness of the presented method for calculating escape zone and the effectiveness of the proposed evasion guidance scheme against multiple pursuing spacecraft.


I. INTRODUCTION
As the number of near-Earth space objects increases, there is growing attention to the importance of the orbital evasion problem.The problem of orbital evasion originates from the problem of orbital rendezvous and interception, in which the spacecraft needs to evade other maneuvering or nonmaneuvering space objects that may be about to collide through its own control [1].
Early studies on orbital evasion were mostly carried out against nonmaneuvering space objects.In these research, the evasion strategies were normally generated by maximizing an optimization index, such as the terminal miss distance [2], [3] or the collision probability [4].However, these evasion strategies solved by the one-sided optimization cannot be used against a maneuvering space object since the possible maneuvers of the space object were not considered in these research.
Taking into account the possible maneuvers of the space objects, the traditional collision avoidance problem will become an orbital pursuit-evasion (OPE) problem.The differential game method [5] with maneuver assumptions for both pursuer and evader is more suitable to deal with this two-sided optimal control problem.
Since the differential game theory [6], [7] was put forward, many studies on OPE games have been carried out.Analytical methods, for example, the closed-form solution of barrier [8], and numerical methods, such as the semidirect collocation nonlinear programming method [9], the multiple shooting method [10], and the dimension-reduction method [11], are two widely used types of methods for solving an OPE game.Moreover, the nonlinear control for OPE game was realized using the state-dependent Riccati equation method in [12].Two optimal guidance methods, namely, the Cartesian model and the spherical model, for the long-distance OPE game were proposed in [13].The PE differential game for satellites with continuous thrust was investigated from the viewpoint of reachable domain in [14].Many other works have focused on solving the saddle point solution of the game quickly and efficiently [15], [16], [17], [18], as well as the more complex dynamic environment [19], [20] and information structure [21].
It should be noted that all of the research abovementioned focused on the two-player OPE games involving one evading spacecraft and one pursuing spacecraft with continuous thrust [22].Almost few studies have been conducted on the orbital evasion strategy against multiple pursuers.Meanwhile, considering the fact that impulsive thrust is still the main form of spacecraft maneuver nowadays, it is required to investigate the impulsive evasion strategy.Chandrakanth and Scheeres [23] conducted research aimed at the two-player impulse OPE game, while the assumption of two-impulse transfers is not suitable for short-distance PE.Overall, there is still a gap in the research concerning the impulsive evasion strategy for the short-distance scenario against multiple pursuing spacecraft.
According to the initial relative distance between the pursuing spacecraft and evading spacecraft, the orbital evasion problems can be divided into two scenarios [23]: long-distance evasion scenario and short-distance evasion scenario.In contrast, the short-distance evasion scenario, which is closely related to the relative motion state between the pursuing spacecraft and evading spacecraft, is the final stage of the complete orbital evasion process, and has the characteristics of strong adversarial and high real-time requirements.This article works on the short-distance orbital evasion problem against multiple pursuers with impulsive thrust.There are two main contributions of this article.The first contribution of this work is to propose the concept of escape zone (EZ) for an impulsive orbital evasion problem against multiple pursuers, and present an effective method for calculating the EZ at any decision moment by introducing the approximate relative reachable domain (RRD) ellipsoid for orbital impulsive maneuver.Second, a two-sided optimization process based on the calculated EZ is performed to generate the optimal pursuit-evasion strategies, which provides a possible idea for solving impulsive orbital games.
The rest of this article is organized as follows.Section II introduces the approximate RRD, and presents the concept of EZ in this article.A method based on the approximate RRD is proposed in Section III to calculate the defined EZ.In Section IV, an escape-zone-based optimal evasion guidance law is proposed based on the escape value (EV) defined through the two-sided optimal strategies analysis.Numerical examples are provided in Section V. Finally, Section VI concludes this article.

A. Relative Reachable Domain Approximate Ellipsoid
The RRD concept is introduced in this article to determine the potential relative state reach set of a spacecraft under given maneuverability.Different from the concepts in [24] and [25], an RRD denoted by D(X 0 , V max , t) here is described as the set of all possible relative positions that can be transferred to after time t from initial relative state X 0 under the maximum available velocity increment V max .Assuming that the reference orbit is circular, the RRD can be established based on the linearized dynamics, i.e., the CW equations [26], which have been proved sufficiently accurate to model the close range relative motion [25], as ⎧ ⎨ ⎩ (1) where and s = sin(ω t), c = cos(ω t), ω is the mean angular motion of the reference orbit, x 0 , y 0 , z 0 , ẋ0 , ẏ0 , and ż0 are components of the initial relative position vector and relative velocity vector in the reference local-vertical, localhorizontal (LVLH) frame, ξ V and η V are two angles characterizing the direction of the initial impulse vector in the reference LVLH frame.
Obviously, the envelope of the RRD described by ( 1) is exactly the boundary of the time-constrained reachable relative state distribution mentioned in [24] and [25] at given time t, without considering the initial position uncertainty.Specifically, the commonality between them is that the fixed-time RRD considered in [24] and [25] will be the same as that in this article if the initial position uncertainties r are set to 0 and the velocity increment applied at the decision moment is considered as the initial velocity uncertainties v.The difference between them is that the work presented in [24] and [25] focuses on solving the accurate inner and outer boundaries of the RRD in any direction within a relative motion time range, while this article focuses on efficiently solving the approximate envelope of the RRD at any given time for quick analysis of the positional relationship between RRDs for different participants in the OPE games.
When the relative motion time is much smaller than the period of the reference orbit, the envelope of the RRD at time t can be approximated as where It should be noted that the approximation in (4) occurs only in the x-y reference plane not in the z direction.Specifically, the right side of the first formula in (1) is approximated by(κ xy V max cos η V ) 2 .Let F 1 ( t) and F 2 ( t) denote the former and latter expression, respectively.F 1 (0) = F 2 (0) and lim t→0 |(F 2 − F 1 )/F 1 | = 0 can be proved, which demonstrate the equivalence of F 1 and F 2 at t = 0. Furthermore, a quantitative analysis is given as follows to investigate that this approximation is valid for how small the motion time is.
Suppose that the height of the reference orbit is 35 786 km, and that V max = 10 m/s.For each time t taken from 0 s to 7% times the reference orbital period T r , a total of 1 000 000 Monte Carlo runs calculating the absolute value of relative error (ARE) of relative reachable The results of the Monte Carlo runs are shown in Fig. 1.As seen in Fig. 1, the mean ARE is less than 3% (considered as the maximum value for an acceptable approximation error) when t is less than about 6.665%T r = 5741.372s.Similarly, the maximum ARE is less than 3% when t is less than about 4.797%T r = 4132.237s.The results for other orbital heights and impulse magnitudes are also analogous to that in Fig. 1, some of which are given in Tables I and II.
Accordingly, the following conclusions are obtained about the approximate envelope defined in (4).
1) The approximate envelope of the RRD at time t is an ellipsoid of revolution (termed RRDE) in the reference LVLH frame, which is centered on the nominal relative position at time t, with a = b = κ xy V max and c = κ z V max as the major and minor semi axes, respectively, as shown in Fig. 2.
2) The error caused by the approximation is acceptable only when t is much smaller than T r (termed shortterm assump-tion).Specifically, t needs to be less than about τ 1 = 6.665%T r or τ 2 = 4.797%T r so as to ensure that the mean ARE or maximum ARE is less than 3%, respectively.3) Under the short-term assumption, a and c both increase monotonically with increasing time t.

B. Problem Statement and Escape Zone Definition
Consider an orbital evasion problem against N pursuing spacecraft.Each participating spacecraft is controlled by a three-dimensional impulse input at each decision moment, which represents the thrust effect during this decision period.The pursuing spacecraft attempts to satisfy some certain intercept conditions, for example, position matching with the evader.Contrarily, the evader needs to avoid entering this terminal interception range through its own control.
Before continuing, the following assumptions are made.
1) All the participating spacecraft are assumed to move in a two-body gravitational field, without considering any forms of perturbations and real-world noises.
2) The maneuverability of the evader is assumed to be stronger than that of the pursuer to ensure escape is possible [27].Specifically, two sides of the game are assumed to maneuver at the same time [28].Meanwhile, the decision periods for the pursuers and evader are the same, but the maximum available velocity increment of the single impulse for the evader is slightly larger than that of the pursuer.
3) The observation information is assumed to be complete, that is, each participating spacecraft can accurately acquire the status of other spacecraft in real time.
As shown in Fig. 3, the evasion process starts from the moment t 0 when the minimum distance R min between the pursuers and evader is less than a given alert distance R alert , for example, 100 km, and terminates when R min is less than the minimum valid interception distance R intercept of pursuers (evasion failure) or when R min is larger than R alert again (evasion success).Therefore, the entire pursuitevasion motion can be considered as the close range relative motion.For the above orbital evasion problem, a unified time RRD ellipsoid of the evader (termed UTRRDE) at time t U is proposed first, where t U is a unified time introduced to unify different RRD ellipsoids for different intersection moments herein.There is no strict constraint on the selection of the unified time t U as long as it is larger than the decision moment t 0 when the impulse is applied, i.e., the relative motion time (t Ut 0 ) is positive, to ensure that the geometric parameters a and c of UTRRDE are larger than 0. However, considering the approximation error caused by the RRD ellipsoid, it is suggested to select a unified time less than the maximum valid relative motion time τ 1 = 6.665%T r determined in Section II-A to reduce the guidance error.
In this article, the corresponding point M TU on the UTRRDE for a point M on the RRDE E at time t is defined as the terminal relative position point after time t U on the impulsive relative trajectory (termed the characteristic trajectory), along which the evader can transfer to M after time t.
The corresponding point M TU 1 of a point M 1 on the intersection curve I 1 of two RRDEs at time t can be calculated as follows: The velocity increment of the characteristic trajectory for M 1 and M TU 1 is computed first by where ) represents the state transition matrix from the initial impulse to the relative position at time t in the reference orbital frame, and x M 1 , y M 1 , and z M 1 are components of the relative position vector of M 1 in the reference LVLH frame.Then, the corresponding point M TU 1 can be obtained by In this regard, the concept of (EZ here proposed for the above-mentioned evasion problem is defined as follows: Let t U t denote the EZ at decision moment t with respect to the unified moment t U , where t U > t.Define the initial reference orbit with the initial orbit of evader.As shown in Fig. 4, t U t is a combined area on the RRDE at time (t U -t), and all transfer trajectories that can reach to this area after time (t U -t) through a single velocity impulse, the magnitude of  which is not larger than V max , can escape the interception of pursuers.

III. ESCAPE ZONE CALCULATION
In this section, a method for finding the intersection of two RRDEs is presented first.The largest threatened area on the RRDE of evader by each pursuing spacecraft is solved then, and the EZ is eventually obtained by successively removing these areas from the entire ellipsoid.

A. Finding Intersection of Two RRDEs
Let the initial time t 0 = 0. Define the initial reference orbit with the initial orbit of evader, in which way the center of the RRDE of evader will always be located at the origin of the reference orbital frame.Suppose that the RRDE of pursuer P (RRDE P ) intersects that of evader E (RRDE E ) at time t, as shown in Fig. 5.The intersection of two RRDEs represented by the red line I 1 in Fig. 5 can be solved as follows.
Let O P and O E denote the center of the RRDE P and RRDE E , respectively.For the plane Q 0 determined by the line O P O E and the Z-axis of the reference frame, the components of its normal vector n 0 in the reference LVLH frame are where x cP , y cP , and z cP are the components of the nominal relative position vector of P at time t in the reference LVLH frame.
The plane Q defined by is obtained by rotating Q 0 about the line O P O E by angle α.
According to Rodrigues' rotation formula [29], the parameters p, q, and r in ( 10) can be computed by (11) Let a P , c P and a E , c E denote the geometric parameters of RRDE P and RRDE E , respectively.Then, the ellipsoids of P and E can be defined by The plane Q intersects the RRDE P and RRDE E at the curve I P and I E , respectively, as shown by the green lines in Fig. 5. Evidently, the intersection curves I P and I E must be ellipses or circles [30].If and only if the plane Q is the equatorial plane of the RRDE, the intersection curve is a circle.Let By using the Householder transformation [31], the parametric equation [32] of the intersection curve I i can be given by where x cE = y cE = z cE = 0, and ϕ i ∈ [0, 2π ] is the phase parameter of each point on I i .Then, two intersection points M 1 and M 2 of I P and I E , both of which lie on the intersection curve I 1 of two RRDEs, can be calculated by solving the equations x IP = x IE , y IP = y IE , and z IP = z IE .However, this problem is difficult to solve analytically because of the high nonlinearity in (14).In this regard, the numerical solution is obtained here by minimizing the objective function that for design variables ϕ P ∈ [0, 2π ] and ϕ E ∈ [0, 2π ] through the nonlinear optimization method, for example, the Quasi-Newton method [33].
In order to improve the convergence of the abovementioned optimization process, the phases of the intersection points M 1 and M 2 of two circles obtained by the intersection of plane Q and each RRD approximate sphere (RRDS) are used to provide initial guesses forϕ i , as shown in Fig. 6, where the RRDS is a sphere whose center coincides with the center of RRDE and whose radius is equal to a of RRDE.Meanwhile, the Jacobian matrix of the objective function is provided by where and × {a 2  P pq cos ϕ P − [λ P (c P r {[λ P (c P r − λ P ) + (a P q) 2 ] cos ϕ P − a 2 P pq sin ϕ P } + c P λ P 2c P λ P (a P p cos ϕ P + a P q sin ϕ P ) + a E q sin ϕ E ) (a P q cos ϕ P − a P p sin ϕ P ).(17) In this manner, the intersection curve I 1 of the RRDE P and RRDE E at time t can be obtained ultimately by calculating the intersection points of I P and I E for every rotation angle α from 0 to π. Considering that the three main axes of RRDE P and those of RRDE E are parallel respectively, the symmetry can be utilized here to calculate the intersection points of I P and I E for α from π/2 to π by directly using the results for α from 0 to π/2 according to ⎧ ⎨ ⎩ where M 4 (x 4 , y 4 , z 4 ) is the symmetrical intersection point for α from π/2 to π of the intersection point M 3 (x 3 , y 3 , z 3 ) for α from 0 to π/2.Therefore, the computational effort can be eventually reduced by only calculating the intersection points of I P and I E for α from 0 to π/2.
Because the calculated intersection curves I 1 for different intersection times are located on different RRDEs, it is hard to compare them quantitatively.To solve this problem, the unified time RRD ellipsoid presented in Section II-B is introduced here.Through ( 6)-( 8), the intersection I 1 for any intersection time can be projected onto the same UTRRDE to obtain the corresponding curve I 2 , as shown in Fig. 7.

B. Solving the Largest Coverage on UTRRDE
To solve the largest coverage on the UTRRDE swept by the pursuer P during the entire relative motion process, the RRDE of P and the RRDE of the evader E are propagated along the time axis.First, the positional relationship between two RRDEs at the given time t can be judged as follows.
Rewrite (12) into the following quadratic equations: where Since A E and A P are positive definite, Q i (X) < 0, (i = P, E) defines the inside of the RRDE i , and Q i (X) > 0 correspond to the outside.
Generally, Q E (X) = λ define all level curve ellipsoids [34] of RRDE E , where the minimum (negative) value of λ and λ = 0 define the O E and the RRDE E , respectively.The minimum level value λ 0 and maximum level value λ 1 of the level curve ellipsoids that intersect the RRDE P can be calculated by minimizing Q E (X) subject to the constraint Q P (X) = 0, which can be solved by the method of Lagrange multipliers [35].
Define the Lagrange function as where ε is the introduced Lagrange multiplier.Differentiating (21) yields where the gradient ∇ represents the derivatives in X.
Setting (22) equal to zero yields Formally solving for X yields where δ(ε) = det(A E + εA P ) is a cubic polynomial in ε, and Y(ε) is a column vector that has components cubic in ε.
Setting (23) equal to zero and replacing X with (25) yields which is a degree six polynomial in ε.The X can be computed by substituting the computed roots [36] of ( 26) into (25), and the corresponding value of Q E (X) can be calculated by (19).The minimum and maximum values of the calculated Q E (X) are exactly λ 0 and λ 1 , respectively.The positional relationship between two RRDEs at time t can be decided by comparing λ 0 and λ 1 with 0: If λ 0 > 0, two RRDEs are separated.If λ 0 < 0 and λ 1 > 0, two RRDEs intersect.If λ 1 < 0, the RRDE P is contained in the RRDE E .
Then, three cases about the types of the nominal relative trajectory of P are discussed according to the positional relationships between two RRDEs along the time axis.Let T P denote the largest coverage of P. As mentioned, for the separation case, T P = ∅.For the intersection case, T P is the union of all the areas enclosed by the intersection curve I 2 during the time interval (t ex1 , t ex2 ), as shown in Fig. 9(a).Similarly, for the contained case, T P is the union  of all the areas enclosed by the intersection curve I 2 during the time interval (t ex1 , t in1 ) and (t in2 , t ex2 ), as shown in Fig. 9(b).
Considering the infinity of numbers of the rotation angle α around the line O P O E and the intersection moment t on the time axis, an ellipsoidal polygon with finite vertices, which all lie on the intersection curve I 2 , is used to replace I 2 , and a time step is selected for discrete calculations along the time axis.Specifically, the larger the number of polygon vertices and the smaller the time step, the higher the solution accuracy but the higher the computational burden.Under this approximation, the union of two intersection curve areas at different intersection time can be calculated by judging the inclusion relationships of the vertices of two ellipsoidal polygons [37].
In this way, T P of each pursuing spacecraft is recorded as a data table, which takes the intersection time t as the data index and the geodetic coordinates of all vertices of the approximate polygon of the intersection curve I 2 at time t as the data value, and is stored in a database for all the pursuers.The boundary of T P is recorded as a series of vertices associated with the intersection times.
Eventually, the EZ can be obtained by considering the threats of all the pursuing spacecraft.Let T Pi , (i = 1, 2, …, N) denote the largest coverage of pursuing spacecraft P i solved at the decision moment t.For R intercept = 0, the EZ t U t (see Fig. 10) can be calculated ultimately by removing all these largest coverage from the entire unified time ellipsoid where 0 is the entire UTRRDE at the unified time t U .

IV. OPTIMAL EVASION GUIDANCE LAW
In this section, the analysis of the two-sided optimal strategies for this orbital evasion problem is presented first, and the concept of the EV is proposed next.Finally, a proposed EV-optimal guidance law is given to avoid interception by multiple pursuing spacecraft.

A. Two-Sided Optimal Strategies Analysis
In order to ensure the security of the evasion trajectory, the optimal interception strategy of the pursuers for a given evasion trajectory needs to be analyzed.
Before discussion, the concepts of the evasion impulse and the pursuit impulse as shown in Fig. 11 are clarified first.An evasion impulse here denoted by V E [r d ,t U ] represents the impulse of the characteristic trajectory with the terminal position r d on the UTRRDE.A pursuit impulse here denoted by V Pi [r d ,t] represents the impulse of the relative transfer trajectory of P i that can make P i intercept the evader after time t, which moves along the characteristic trajectory with the terminal position r d on the UTRRDE.According to Section III, the terminal position r d is exactly located on the intersection curve I 2 of the intersection time t on the UTRRDE.
Two cases about the two-sided optimal strategies at decision moment t are discussed as follows.
1) If t U t = ∅, i.e., the calculated EZ is not empty as shown in Fig. 12, then any evasion impulse V E [r d ,t U ] satisfying that r d ∈ t U t leads to a successful evasion, which means that the pursuers cannot intercept the evading spacecraft with abovementioned evasion impulse in a limited time, no matter what strategy they adopt.
In this case, the pursuers will attempt to minimize the miss distances with the maneuvered evader intuitively.On the contrary, the evading spacecraft needs to maximize the minimum miss distance of all pursuers among all the successful evasion trajectories.
The miss distance M Pi is defined as the minimum distance between the maneuvered pursuer P i and the maneuvered evader E along the time axis, which can be calculated by M Pi = min{|r Pi (t ) − r E (t )||t ∈ [0, τ ]}, where r Pi (t) and r E (t) represent the relative position vector of the maneuvered pursuer P i and the maneuvered evader E at time t in the reference LVLH frame, respectively.Obviously, the optimal strategy of P i is to adopt the maneuver impulse that can minimize M Pi .Actually, according to the definition of the largest coverage on the UTRRDE, the minimum miss distance must fall on the situation in which the terminal position of the characteristic trajectory is located on the largest coverage T Pi since T Pi stands for the boundary of the maximum available velocity increment of P i .Therefore, the optimal pursuit impulse V Pi [r d2 ,t] can be solved by searching r d2 and t on T Pi to minimize the miss distance M Pi , which can be represented as an optimization problem as follows: min (28) The minimum miss distance of all pursuers is finally obtained by M min = min{M P1 , M P2 , . . ., M PN }.As mentioned, the optimal evasion impulse V E [r d1 ,t U ] can be obtained by searching r d1 on the EZ to maximize M min , which can be represented as min 2) If t U t = ∅, i.e., the EZ does not exist as shown in Fig. 13 because the UTRRDE has been completely covered by the threat areas of all pursuers, then the pursuers can always find an appropriate strategy to intercept the evader, no matter what maneuver impulse the evader adopts.
In this case, the pursuers will attempt to minimize the interception time for each given evasion impulse while the evading spacecraft attempts to maximize the minimum time to be intercepted among all the possible evasion trajectories.
With the help of the established database in Section III, the minimum time t min Pi to be intercepted by P i for a given evasion trajectory pointing to a certain terminal position r d2 on the UTRRDE can be computed by searching in the database for the minimum intersection time satisfying that the terminal position r d2 is inside the area enclosed by intersection curve I 2 at this intersection time.Then, the optimal impulse of P i can be determined as V Pi [r d 2 , t min  Pi ].Let t min P denote the minimum time among all the valid t min Pi .The optimal evasion impulse V E [r d1 ,t U ] can be obtained by searching r d1 on 0 to maximize t min P , which can be represented as min B. Escape-Value-Optimal Guidance Law According to the two-sided optimal strategies, the EV is defined here as a scalar describing the effect to escape of an evasion trajectory pointing to the terminal position r d on the UTRRDE, which can be computed by where t min P is the minimum time to be intercepted of the evader for the given evasion impulse V E [r d ,t U ] when the EZ is empty, T SI is the dimension of time, M min is the minimum miss distance of all pursuers for the given evasion impulse V E [r d ,t U ] when the EZ is not empty, and L SI is the dimension of distance.
The EV described previously can be used as the basis for a closed-loop evasion guidance scheme for an evading spacecraft against multiple pursuing spacecraft.
In this scheme, at each decision moment of the evader, the evasion impulse is determined to obtain the largest EV, which generally implies the greatest evasion possibility.Corresponding to two different cases in the EV calculation, when the EZ exists, the optimal evasion impulse is desired to achieve the largest minimum miss distance of all pursuers.In addition, when the EZ does not exist, the evasion impulse is still optimized for acquiring the maximum time to be intercepted by the pursuers, which means that every nonoptimal maneuver of the pursuing spacecraft will result in the increase of the interception time and may eventually lead to successful evasion.Such a guidance scheme could be implemented as follows.
1) Obtain the current relative states of all N pursuers.2) For the pursuing spacecraft P i , propagate the nominal relative motion from the current state to time T r /2 using the nonlinear relative dynamic model [38], which is a fairly sufficient choice to ensure that all externally and internally tangent moments of two RRDEs can be taken into account.3) Choose a small time step t s for the RRDE propagation.Calculate λ 0 , λ 1 and judge the positional relationship between the RRDE Pi and RRDE E along the nominal trajectory by ( 19)- (26).Determine the type (separation, intersection, or contained) of the nominal trajectory for P i as in Fig. 8, and record the time t ex1 and t ex2 for the intersection case or time t ex1 , t in1 , t in2 and t ex2 for the contained case.
4) For the separation case, let the largest coverage T Pi be empty, and skip to step 7).5) For the intersection or contained case, calculate all the intersection curves I 1 of the RRDE Pi and RRDE E from 0 to τ 1 with time step t s using the method in Section III-A, and project them onto the UTRRDE through ( 6)-( 8) to obtain the corresponding curves I 2 .Save all the curves I 2 into the database for P i .6) Calculate the union of all I 2 to obtain T Pi .7) Repeat the procedure from steps 2) to 6) until the largest coverage of all pursuers has been calculated.8) Calculate the EZ through (27).9) Determine the optimal evasion maneuver by searching r d on the EZ (if the EZ is not empty) or on 0 (if the EZ is empty) to minimize the EV computed by using (31).10) This procedure, given by the preceding steps, is continued until the terminal conditions of the evasion process are satisfied.If R min is larger than R alert , then successful evasion has been attained.
It should be noted that this guidance scheme provides a conservative evasion strategy for the evading spacecraft, in which all the pursuing spacecraft are considered "smart" enough, i.e., will adopt the optimal pursuit strategy.In other words, the provided evasion strategy is a greedy and locally optimal solution rather than the global optimal solution for the entire game with more than one round.However, thanks to the calculated EZ, this local optimal evasion strategy can be quite concise, fast, and effective.In addition, considering the uncertain antagonism of the pursuers, it is extremely difficult to accurately determine the total number of rounds of the entire game when the strategies of pursuers are unknown, while the guidance scheme could guarantee the basic proceeds of the evading spacecraft in the real-time games.

V. NUMERICAL EXAMPLES
In this section, the effectiveness of the proposed evasion guidance scheme is verified by a numerical simulation of an orbital evasion scenario against several pursuing spacecraft.
An evading spacecraft E and four pursuing spacecraft P 1 , P 2 , P 3 , P 4 are included in this orbital evasion scenario, which starts from the initial time t 0 = 0 s when the minimum distance R min between the pursuers and evader is less than the given alert distance R alert = 100 km of evader.The orbital elements at t 0 , including semimajor axis a 0 , eccentricity e 0 , inclination i 0 , right ascension of ascending node 0 , argument of perigee ω 0 , and mean anomaly M 0 , are given in Table III.
The parameters of mass and maneuverability of the spacecraft are given in Table IV.As mentioned in Section II-B, the minimum time intervals T d between two impulses of pursuers and that of the evader are set to be equal, and the maximum available velocity increment 2.4 m/s for each impulse of the evader is set to be larger than that 2 m/s of the pursuer, which conforms to the   assumption presented in Section II-B that the maximum available velocity increment of the single impulse for the evader is slightly larger than that of the pursuer.
Define the initial reference orbit with the initial orbit of evader.The initial nominal trajectories of pursuers in the reference LVLH frame are portrayed in Fig. 14    to Table III.Suppose that the minimum valid interception distance of pursuers R intercept = 1 km.Then, the initial interception time along the nominal trajectory can be computed as about 4296.284s for P 3 and 4295.958s for P 4 , which indicates the necessity of evasion for E.
Two cases of different pursuing strategies are considered: case 1) the optimal pursuing strategy described in Section IV-A; and case 2) the minimum zero effort miss interception strategy solved by the differential evolution algorithm [39] at each decision moment with the constraint of the maximum available velocity increment V maxP for each impulse.
The evading spacecraft is guided through the proposed guidance scheme in both cases.All the spacecraft are simulated to move in a two-body gravitational field in the numerical examples.The results and analysis of the numerical examples are presented as follows.
For case 1, the evasion trajectory of evader and the pursuit trajectories of pursuers in reference frame are portrayed in Fig. 15    portrays the time history of the distance between each pursuing spacecraft and the evading spacecraft.As shown in the figures, the evading spacecraft guided through the proposed scheme successfully escape from all the pursuing spacecraft at 15715.374 s, at which R min is larger than R alert again.
It is shown in these figures that the entire evasion process can be approximately considered as two stages separated at about 6000 s in case 1.In the first stage, the evading spacecraft attempted to break through the encirclement of the pursuing spacecraft, after which the evader focused on keeping away from all the pursuers through its stronger maneuverability in the second stage.The relative position history shown in Fig. 20 and the pursuit-evasion velocity increment history in the inertial frame shown in Fig. 21 also confirm this ratiocination.The detailed transfer trajectories within 6000 s are portrayed in Fig. 22, where the dashed lines denote the initial nominal trajectories of the pursuers.
All the intersections I 2 during the intersection time intervals and the largest coverage T P1 calculated for P 1 at t 0 in the reference LVLH frame are shown in Fig. 23    all I 2 and T P1 on the UTRRDE.Similarly, the intersections and the largest coverage for P 2 , P 3 , and P 4 at t 0 are portrayed in Figs.24-26, respectively.As shown in these figures, the results of the largest coverage coincide with the conclusion obtained in Section III-B.As a result, the EZ at t 0 is ultimately obtained through (27).The geodetic coordinates of the EZ at t 0 are shown by the gray area in Fig. 27.Obviously, the EZ is equal to the remaining area after removing all the calculated largest coverage of pursuers from the entire UTRRDE.In the above simulations, the number of polygon vertices to replace I 2 was set to 360 with a step of 1°for the rotation angle α, and the time step for discrete calculations was set to approximately 5.562 s corresponding to the intersection time interval equally divided into 200 segments, which is preliminarily determined by the intersection time range of two RRDSs.Under such parameter configurations, the computing time for a data table is about 10.934 s on an Intel(R) Core(TM) i7-7700 on a 3.60 GHz CPU.
As for the two-sided optimal strategies, considering that the EZ exists at t 0 , the optimal evasion impulse is desired to achieve the largest minimum miss distance M min of all pursuers.According to Section IV-A, the maximum M min is obtained as 1.596692 km when the longitude of the terminal evasion position r d on UTRRDE is −3.376190°and the latitude is 37.528571°.Correspondingly, Fig. 28     = −3°and latitude = 38°, which is close to the optimization results.
For case 2, the evasion trajectory of evader and the pursuit trajectories of pursuers are portrayed in Figs.29-32.Fig. 33 portrays the time history of the distance between each pursuing spacecraft and the evading spacecraft.As shown in the figures, the evading spacecraft successfully escape from all the pursuing spacecraft at 14353.583 s, which is smaller than 15715.374s in case 1.In addition, as seen in Fig. 33, the minimum distance between pursuers and   evader during the entire evasion process is 13.020359 km in case 2, which is larger than that of 6.462936 km in case 1.This is intuitive and can be accounted for by the nonoptimality of the pursuing strategy in case 2 since the pursuers did not consider any possible maneuver of the evader in this case.
It should be noted that although the total time of the evasion process is relatively larger than τ , the effect of guidance can still be guaranteed because the reference orbit for optimal impulse calculation is updated to the current orbit of evader at every decision moment to make sure the RRDE of evader is always located at the origin of reference frame while calculating the EZ, and the maneuver period T d is much smaller than τ .

VI. CONCLUSION
In this article, an escape-zone-based optimal evasion guidance law was proposed for an evading spacecraft on near-circular reference orbit against multiple pursuing spacecraft.The concept of EZ was presented for the impulse orbital evasion problem first.A method based on the approximate RRD ellipsoid was proposed to calculate the EZ against multiple pursuers for each maneuver moment.The EV was defined according to the two-sided optimal strategies analysis and eventually used as the basis for the closed-loop evasion guidance scheme.
The proposed evasion guidance scheme provides a conservative EV-optimal evasion strategy for the evading spacecraft by considering the optimal intercept strategies of pursuers such that to obtain the largest minimum miss distance of all pursuers when EZ exists or the maximum time to be intercepted by pursuers when EZ is empty.Finally, the numerical examples of an orbital evasion scenario with one evader and four pursuers were provided to verify the effectiveness of the presented method for EZ calculation and the proposed guidance scheme.The simulation results showed that the evader successfully escaped from all the pursuers in both cases against the optimal and nonoptimal pursuing strategies.Moreover, a better performance of evasion process, including a larger minimum distance and a smaller evasion time, against the nonoptimal pursuing strategy was found than that against the optimal pursuing strategy.

APPENDIX A. Derivation for Approximate RRD
Considering a short-distance relative motion scenario with the circular reference orbit, after time t of applying the impulse V with the components [ V x V y V z ] in the reference orbital frame, the components of the relative position of the spacecraft based on CW equations satisfy where ω is the mean angular motion of the reference orbit, x 0 , y 0 , z 0 , ẋ0 , ẏ0 , and ż0 are components of the initial relative position vector and relative velocity vector in the reference LVLH frame, and ξ V and η V are two angles characterizing the direction of the initial impulse vector in the reference LVLH frame.
Move the terms related to the initial state in ( 32) and (33) to the left side of the equations.Taking the square of both sides and adding the two equations of x and y gives Then, (34) and (35) can respectively be written as Equations ( 38) and (39), i.e., (1) in Section II-A, are the calculation formulas for the 3-D RRD at time t.Let F 1 (t) denote the right side of (38).When the relative motion time t is much smaller than the period T r of the reference orbit, use F 2 (t ) = [κ 2 + (κ 1 + κ 3 + κ 4 + κ 5 )/2]( V cos η V ) 2 to approximate F 1 (t).In this way, the approximate envelope of the RRD can be obtained as The rationality of the above approximate process is proved from the following two aspects.First, the equivalence of F 1 (t) and F 2 (t) in the neighborhood of t = 0 for the above-mentioned approximation process is demonstrated as follows.

B. Guidance Scheme for Pursuing Spacecraft in Simulation Case 1 of Section V
The guidance scheme for the pursuing spacecraft P i in simulation case 1 of Section V could be implemented as follows.
1) Obtain the current relative state of the evading spacecraft and calculate the state relative to evader.2) Propagate the nominal relative motion from the current state to time T r /2 using the nonlinear relative dynamic model.3) Choose a small time step t s for the RRDE propagation.Calculate λ 0 , λ 1 and judge the positional relationship between the RRDE Pi and RRDE E along the nominal trajectory by ( 19)- (26).Determine the type (separation, intersection, or contained) of the nominal trajectory for P i as in Fig. 8, and record the time t ex1 and t ex2 for the intersection case or time t ex1 , t in1 , t in2 , and t ex2 for the contained case.4) For the separation case, let the largest coverage T Pi be empty, and skip to step 7).5) For the intersection or contained case, calculate all the intersection curves I 1 of the RRDE Pi and RRDE E from 0 to τ 1 with time step t s using the method given in Section III-A, and project them onto the UTRRDE through ( 6)-( 8) to obtain the corresponding curves I 2 .Save all the curves I 2 into the database for P i .6) Calculate the union of all I 2 to obtain T Pi .7) Determine the optimal pursuit maneuver by searching r d2 and t on T Pi to minimize the miss distance M Pi or directly as V Pi [r d2 ,t min Pi ]. 8) This procedure, given by the preceding steps, is continued until the terminal conditions of the orbital game are satisfied.If R min is smaller than R intercept , then successful interception has been attained.
. Projections of these transfer trajectories on the X-Y plane, X-Z plane, and Y-Z plane of the reference LVLH frame are shown in Figs.16 -18, respectively.Fig. 19

Fig. 19 .
Fig. 19.Distance history between each P and E in case 1.
(a) by the blue dotted lines and the red solid line, respectively.Meanwhile, Fig. 23(b) portrays the geodetic coordinates of

Fig. 23 .
Fig. 23.Calculated I 2 and T P1 for P 1 of intersection type at t 0 .

Fig. 24 .
Fig. 24.Calculated I 2 and T P2 for P 2 of intersection type at t 0 .

Fig. 25 .
Fig. 25.Calculated I 2 and T P3 for P 3 of contained type at t 0 .
illustrates the heat map of M min calculated for all the terminal evasion positions generated by traversing the longitude ∈[−180°, 180°] and latitude ∈[−90°, 90°] with the step of 1°.The maximum M min calculated is 1.576265 km when longitude

Fig. 28 .
Fig. 28.Heat map of M min calculated for different r d on EZ.

Fig. 31 X
Fig. 31 X-Z projection of transfer trajectories in case 2.

Fig. 33 .
Fig. 33.Distance history between each P and E in case 2.

TABLE I
Results for More Reference Orbital Heights

TABLE II Results
for More Impulse Magnitudes distances are performed with different ξ

1 )
Case 1, separation: If λ 0 > 0 holds for any time t, then the largest coverage does not exist since the two RRDEs do not intersect at any time [see Fig. 8(a)].2) Case 2, intersection: If λ 1 > 0 holds for any time t, and there exist time t such that λ 0 < 0, then the two RRDEs intersect in a period of time between t ex1 and t ex2 , where t ex1 and t ex2 represent the first and the second externally tangent time of two RRDEs, respectively [see Fig. 8(b)].Therefore, the largest coverage is a continuous area on the UTRRDE swept by P during the time interval (t ex1 , t ex2 ). 3) Case 3, contained: If there exist time t such that λ 1 < 0, then the two RRDEs intersect in two time intervals between which the RRDE P is totally contained in the RRDE E , as shown in Fig. 8(c).Therefore, the largest coverage contains two separate areas on the UTRRDE successively swept by P during the time interval (t ex1 , t in1 ) and (t in2 , t ex2 ).

TABLE III Initial
Orbital Elements of Spacecraft

TABLE IV Mass
and Maneuverability Parameters of Spacecraft