# Research on Minimum Time Interception Problem with a Tangent Impulse under Relative Motion Models

Authors: Wenzhe Ding; Xinhong Li; Hong Yang
Journal: International Journal of Aerospace Engineering
 (2019)

---

## Abstract

The minimum time interception problem with a tangent impulse whose direction is the same as the satellite’s velocity direction is studied based on the relative motion equations of elliptical orbits by the combination of analytical, numerical, and optimization methods. Firstly, the feasible domain of the true anomaly of the target under the fixed impulse point is given, and the interception solution is transformed into a univariate function only with respect to the target true anomaly by using the relative motion equation. On the basis of the above, the numerical solution of the function is obtained by the combination of incremental search and the false position method. Secondly, considering the initial drift when the impulse point is freely selected, the genetic algorithm-sequential quadratic programming (GA-SQP) combination optimization method is used to obtain the minimum time interception solution under the tangent impulse in a target motion cycle. Thirdly, under the high-precision orbit prediction (HPOP) model, the Nelder-Mead simplex method is used to optimize the impulse velocity and transfer time to obtain the accurate interception solution. Lastly, the effectiveness of the proposed method is verified by simulation examples.


---

## Body

## 1. Introduction

For the orbital interception problem under the two-body model, the Lambert method can be used to solve it when the initial orbital elements of a satellite and a target are known [1]. When the impulse point and the interception point are given, the orbital transfer time can be obtained by the Kepler equation, and the initial velocity required for the orbital transfer can be solved by expressing the transfer time as a univariate function of other parameters [2, 3]. When the relative distance between the satellite and the target is small, the state transition matrix can be constructed to solve the initial velocity required for the orbital transfer. On the relative motion problem, if the target is running on a circular orbit, the Clohessy-Wiltshire (CW) equation can be used to describe the relative motion [4]. For the case where the target is running in an elliptical orbit, the Tschauner-Hempel (TH) equation can be used to describe the relative motion [5]. In 2002, Yamanaka and Ankersen obtained the state transition matrix described by the true anomaly by solving the TH equation [6]. Therefore, for the orbital interception problem in relative motion, the state transition matrix can be utilized to solve the impulse velocity required by the interceptor at the impulse moment.For orbital interception tasks such as space debris removal, the tangent impulse has a simpler impulse direction, which makes the attitude adjustment of the satellite at the impulse point more convenient. So tangent impulse interception is a better interception method. For the tangent problem of coplanar elliptical orbits, Adamyan et al. solved the cotangent transfer orbit by the geometric method and obtained the analytical expression between the orbital parameters and the velocity vector [7]. Thompson et al. studied three types of tangent problems by using the Hodograph theory [8]. However, to ensure that the interceptor and the target have the same flight time, the above research does not apply to the problem of orbital interception. In 2012, Zhang et al. obtained the conditions for the existence of transfer solutions for three types of tangent orbits by using the relationship between the orbital semilatus rectum and the flight direction angle [9]. At the same time, the tangent impulse intercept problem with the minimum time when the target orbit is an ellipse is also studied [10]. Wang et al. studied this when the orbit was hyperbolic [11]. However, the above research contents are all based on absolute motion. For the unguided close-range interception problem requiring shorter interception time, it is necessary to study the tangent impulse interception with the minimum time under the condition of relative motion.In addition to the above studies, there are many related studies on the tangent impulse orbit maneuver problem in recent years [12–15]. But for these studies, the models used are simplified models. Therefore, in order to reduce the impact of the perturbation on the intercepting orbit, it is necessary to further optimize the interception orbit. For the Lambert problem considering J2 perturbation, under the premise of setting the terminal precision, the state transition-sensitive matrix of the two-body model is often used to iteratively obtain the required initial velocity by the shooting method [16]. When the number of flight laps is vast, the homotopy method can be used to divide the entire time interval into small intervals so that the initial velocity can gradually converge [17]. However, for the high-precision extrapolation model adopted in this paper, the above shooting method is no longer applicable due to the lack of useful gradient information. At the same time, tangent interception limits the direction of impulse velocity, which makes it impossible to solve the interception problem by fixing the terminal position.The minimum time interception problem can be classified as an optimization problem. Therefore, the corresponding optimization model can be established by combining different optimization indicators and solved by a direct method or indirect method. The direct method transforms the optimization problem of intercepting orbit into a nonlinear programming (NLP) problem and uses the optimization algorithm to solve it. In recent years, some scholars have used hybrid optimization algorithms to solve the single-impulse interception orbit optimization problem [18, 19]. Among them, GA, as a global optimization algorithm, is insensitive to initial values and has strong robustness, which can exhibit strong global search ability. The disadvantage is that the local search ability is weak and the result precision is low. Therefore, GA often provides initial values for gradient information-based optimization methods [20]. The SQP method is sensitive to the initial value, has a small convergence radius, and is easy to fall into the local optimum. But for the NLP problem, it can quickly converge to get a high-precision solution [21]. Therefore, combining GA and SQP is an effective method for solving orbit optimization problems.Based on the above analysis, the paper studies the minimum time interception problem when the tangent impulse is used. First of all, the relative motion model based on an elliptic orbit is used to transform the interception solution into a univariate function only about the true anomaly of the target. Next, all the solutions in the feasible region are obtained by numerical iteration. Then, considering the initial drift segment, the GA-SQP combination optimization method is used to obtain the minimum time interception solution under the tangent impulse in a target motion cycle. Finally, in the absence of effective gradient information, the Nelder-Mead simplex method [22] is adopted to optimize the impulse velocity increment and the transfer time to obtain an accurate interception solution under the high-precision extrapolation model. The main innovation of this paper is to provide an accurate method to solve the interception solution of the minimum time tangent impulse under relative motion by combining analytical, numerical, and optimization methods. The advantages of this method are as follows: (1) Combine the GA and SQP method to optimize the impulse position, and then obtain the accurate minimum time tangent impulse interception solution without providing the initial value. (2) The Nelder-Mead simplex method adopted in the optimization process avoids the dependence on the gradient information. Under the condition that the target is relatively close to the satellite, the accurate solution under the high-precision HPOP model can be obtained just by the initial solution provided by the linear TH equation.


## 2. Existence Condition of Tangent Impulse Interception Solution

For the tangent impulse interception problem of space targets, it is necessary to judge the existence of the interception solution of the satellite firstly at the impulse moment and give the feasible domain of the solution. In this regard, a detailed proof has been given in the literature [9, 10]. According to the conclusion in the literature [10], for the shooting point P1 in Figure 1, there is a solution only when the interception point is located at S1P2S2⌢.Figure 1
The existence condition of the tangent impulse interception solution.In the above figure, the orbits of the target and the satellite are both anticlockwise.P1 is the shooting point of the satellite; v is the velocity of the satellite for this moment; S1 and S2 are the two intersections of the satellite’s velocity direction and the target orbit, respectively; P2 is the intercepting point on the target orbit; F1 is the focus of the target orbit as well as the focus of the satellite orbit; RP1 and RS1 are radius vectors of the P1 and S1 points, respectively; γ¯ is the flight-direction angle of the satellite at this moment; γ is the flight-path angle of the satellite at the moment; θ is the transfer angle of the projectile.In order to describe the existence range of the interception solution better, according to the conclusion in the literature [9], the interception solution is expressed as a form of the true anomaly of the target, as shown in
(1)c1cosf2+c2sinf2>c3,where
(2)c1=sinϖ2−ϖ1−f1−γ¯+e2R1p2sinγ¯,c2=cosϖ2−ϖ1−f1−γ¯,c3=−R1p2sinγ¯,where ϖ1 and ϖ2 are the arguments of perigee of the satellite and the target orbit, respectively, f1 and f2 are the true anomalies of the satellite and the target orbit, respectively, e2 and p2 are the eccentricity ratio and semilatus rectum of the target orbit, respectively, and R1 is the distance between the satellite and the center of the earth.
(1)
Ifc3/c12+c22∈−∞,−1, then f2∈0,2π
(2)
Ifc3/c12+c22∈1,∞, there is no solution for f2
(3)
Ifc3/c12+c22∈−1,1, then f2∈f2,s1,f2,s2where
(3)f2,s1=arcsinc3c12+c22−atan2c1,c2,f2,s2=π−arcsinc3c12+c22−atan2c1,c2.When the transfer orbit of the projectile is an ellipse, the constraint is [9]
(4)c4cosθ+c5sinθ>c6,where
(5)c4=p2R1cos2γ¯+e21−cos2γcosϖ2−ϖ1−f1,c5=p2R1sin2γ¯+e21−cos2γsinϖ2−ϖ1−f1,c6=p2R1+cos2γ¯−1,where θ is the transfer angle of the projectile.As the analysis of equation (1), there are three situations to discuss:
(1)
Ifc6/c42+c52∈−∞,−1, then f2∈0,2π
(2)
Ifc6/c42+c52∈1,∞, there is no solution for f2
(3)
Ifc6/c42+c52∈−1,1, then f2∈f2,s3,f2,s4where
(6)f2,s3=arcsinc6c42+c52−atan2c4,c5−ϖ2+ϖ1+f1,f2,s4=π−arcsinc6c42+c52−atan2c4,c5−ϖ2+ϖ1+f1.It can be seen from the above conclusion that the feasible domain of the tangent interception solution corresponding to the satellite impulse point isf2∈f2,s1,f2,s2, in which, when f2∈f2,s4,f2,s2, the flight transfer angle of the projectile is greater than π and the interception orbit is hyperbolic. It can be seen from the literature [2] that when θ>π, the transfer time of the intercepting orbit is less than 0. So for any intercepting orbit, the intercepting solution exists in the range f2∈f2,s1,f2,s4.


## 3. Tangent Impulse Interception Solution of Elliptical Orbit in the Relative Motion Model

### 3.1. Relative Motion Model and State Solution under the Elliptical Orbit

The J2000 coordinate system is commonly used as the inertial coordinate system, and itsX-axis is defined to point to the mean equinox at 2000:01:01:12:00:00. The Vehicle Velocity Local Horizontal (VVLH) coordinate system is commonly used as an orbital coordinate system. Its origin is located at the center of the mass of the spacecraft, the X-axis points to the direction of motion of the spacecraft, the Z-axis points to the center of the earth, and the Y-axis meets the right-hand rule.In the VVLH coordinate system, the relative motion of the spacecraft with respect to the target is [6]
(7)ṙtv̇t=A11tA12tA21tA22trtvt+B1tB2taxayaz,where rt=x,y,zT and vt=ẋ,ẏ,żT; ax,ay,azT is the external acceleration, and the coefficient matrix is
(8)A11t=03×3,A12t=I3×3,A21t=ω22−μR230ω̇20−μR230−ω̇20ω22+2μR23,A22t=002ω2000−2ω200,B1t=03×3,B2t=I3×3,where μ is the gravitational constant; ω2=a2μ1−e22/R22 is the orbital angular velocity of the target, and a2 is the semimajor axis of the target orbit.For impulse interception, the external acceleration is 0. Letf20 be the true anomaly of the target’s initial position, and the state equation can be solved as follows:
(9)Xtf2=Φtf2,tf20Xtf20,where tf20 is the initial moment, tf2 is the terminal moment, Xtf20 is the relative state of the satellite with respect to the target at the initial moment, Φtf2,tf20 is the state transition matrix, and Xtf2 is the relative state of the satellite with respect to the target at the terminal moment. The division is performed in and out of the orbital plane, and the expression of Φtf2,tf20 is as follows [6]:
(10)Φtf2,tf20=Φ11tf2,tf20Φ12tf2,tf20Φ21tf2,tf20Φ22tf2,tf20=O11tf2,tf200O12tf2,tf2000g1tf2,tf200g2tf2,tf20O21tf2,tf200O22tf2,tf2000g3tf2,tf200g4tf2,tf20.The state transition equation in the orbital plane is(11)where
(12)where k=μ/h3/2, h is the orbital angular momentum of the target, ρf2=1+e2cosf2, ρf20=1+e2cosf20, sf2=ρf2sinf2, sf20=ρf20sinf20, cf2=ρf2cosf2, cf20=ρf20cosf20, sf2′=cosf2+e2cos2f2, cf2′=−sinf2−e2sin2f2, and J=k2tf2−tf20=∫f20f21/ρ2τdτ.The state transfer equation outside the orbital plane is(13)where
(14)g1tf2,tf20g2tf2,tf20g3tf2,tf20g4tf2,tf20=ρf20−e2sinf21k2ρf2−1cosf2−f20sinf2−f20−sinf2−f20cosf2−f20ρf200−e2sinf201k2ρf20.It can be seen from equations (11) and (13) that the state transition matrix of the relative motion between the satellite and the target is converted into a form related only to the initial true anomaly f20 of the target and the terminal true anomaly f2.


### 3.2. Tangent Impulse Interception Solution under the Elliptical Orbit

(15)T=−sinucosΩ−cosucosisinΩ−sinusinΩ+cosucosicosΩcosusini−sinisinΩsinicosΩ−cosi−cosucosΩ+sinucosisinΩ−cosusinΩ−sinucosicosΩ−sinusini.First, the relative position and relative velocity of the satellite with respect to the target in the VVLH coordinate system are calculated. The coordinate transformation matrix from the J2000 coordinate system to VVLH is shown in formula (15), where Ω is the right ascension of ascending node, u is the argument of latitude, and i is the orbit inclination.The coordinate system conversion matrix of the satellite and the target calculated by equation (15) is T1 and T2, respectively. Since the target VVLH coordinate system is a rotating coordinate system, the relative state of the satellite with respect to the target can be expressed in this coordinate system as
(16)r=T2R1−R2,(17)v=T2V1−V2−ω2×r,where R1,V1 and R2,V2 are the positions and velocity vectors of the satellite and the target in the J2000 coordinate system, respectively, and ω2× is the skew-symmetric matrix represented by the angular velocity vector of the target VVLH coordinate system relative to the J2000 coordinate system.
(18)ω2×=00−ω2000ω200.Equation (17) can be further expressed as
(19)v=T2V1−V2+ω2z0−ω2x.Therefore, the relative positionrf20 and the relative velocity vf20 at the initial moment tf20 can be obtained according to equations (16) and (19).Without drifting, the satellite shoots at the initial moment, according to equations (9) and (10); the relative state of the projectile with respect to the target at the terminal moment is
(20)rd,f2vd,f2=Φ11tf2,tf20Φ12tf2,tf20Φ21tf2,tf20Φ22tf2,tf20rd,f20vd,f20,where rd,f2 and vd,f2 are the relative position and velocity vectors of the projectile with respect to the target at the terminal moment, respectively, and the constraint of rd,f2=0,0,0T is satisfied for the interception problem; rd,f20 and vd,f20 are the relative position and velocity vectors of the projectile with respect to the target at the initial shooting moment, and rd,f20=rf20 is satisfied. Therefore, the velocity vector of the projectile relative to the target at the initial moment can be obtained as
(21)vd,f20=−Φ12−1Φ11rf20.Therefore, in the target VVLH coordinate system, the impulse vector launched by the satellite is(22)Δv=vd,f20−vf20.According to formula (19), the velocity vector of the projectile at the moment of shooting in the J2000 coordinate system is
(23)Vd,f20=T2,f20−1vd,f20−ω2,f20zd,f200−ω2,f20xd,f20+V2,f20.Then, the satellite’s launch impulse vector can be expressed in the satellite VVLH coordinate system as(24)Δv1=T1,f20Vd,f20−V1,f20=T1,f20T2,f20−1Δv.And its unit direction vector can be expressed as(25)Δv^1=Δv1Δv12=Δv^1xΔv^1yΔv^1z.Since the orbits of the satellite and the target are coplanar,Δv^1y=0 is satisfied.The flight-direction angle of the initial moment of the satellite is expressed as a univariate function form with respect to the true anomaly of the satellite:(26)γ¯f20=π2−arctane1sinf101+e1cosf10,where f10 is the true anomaly of the satellite at the initial moment.In order to get the true anomaly of the target corresponding to the interception point, define the functionF:
(27)F=−Δv^1x−Δv^1ztanγ¯f20.In the formula, when the initial moment is determined, the flight direction angleγ¯f20 is constant. At this time, the function F is a univariate function about the true anomaly f2 of the target.When the satellite shoots in a tangent mode, the shooting impulse direction of the satellite is consistent with the velocity direction of the satellite at the moment, that is,(28)Ff2=0.Then, according to the change of the true anomaly of the target, the total interception time of the projectile from the initial launch moment to the encounter moment is obtained:(29)ttotal=tf2−tf20=1k2∫f20f21ρ2τdτ.For equation (29), the Lobatto quadrature method can be used for the numerical integration calculation [23].


### 3.3. Solving the Tangent Impulse Intercept Solution

Since using equation (28) is difficult in obtaining an analytical x of the equation is obtained by the combination of incremental search and the false position method [23]. The specific steps are as follows:Step 1.
According to the state of the initial impulse of the satellite, the feasible domain of the interception solution obtained from the second section isf2∈f2,s1,f2,s4.Step 2.
Using the incremental search method, divide the intervalf2,s1,f2,s4 into n adjacent cells, f2,s1,f2,s1,1, f2,s1,1,f2,s1,2,…, f2,s1,n‐1,f2,s4, and determine whether the function symbol of each pair of interval endpoints changes.Step 3.
Use the false position method to obtain an accurate interception solution in the interval where the symbol changes.IfFf2,s1,mFf2,s1,m+1<0, let
(30)f2,r=f2,s1,m+1−Ff2,s1,m+1f2,s1,m−f2,s1,m+1Ff2,s1,m−Ff2,s1,m+1.Takingf2,r as the new endpoint of the interval, it is iterated by equation (30). When f2,r−f2,r+1<δ, the iteration ends, and the exact numerical solution of the interception solution in interval f2,s1,m,f2,s1,m+1 is obtained.


## 3.1. Relative Motion Model and State Solution under the Elliptical Orbit

The J2000 coordinate system is commonly used as the inertial coordinate system, and itsX-axis is defined to point to the mean equinox at 2000:01:01:12:00:00. The Vehicle Velocity Local Horizontal (VVLH) coordinate system is commonly used as an orbital coordinate system. Its origin is located at the center of the mass of the spacecraft, the X-axis points to the direction of motion of the spacecraft, the Z-axis points to the center of the earth, and the Y-axis meets the right-hand rule.In the VVLH coordinate system, the relative motion of the spacecraft with respect to the target is [6]
(7)ṙtv̇t=A11tA12tA21tA22trtvt+B1tB2taxayaz,where rt=x,y,zT and vt=ẋ,ẏ,żT; ax,ay,azT is the external acceleration, and the coefficient matrix is
(8)A11t=03×3,A12t=I3×3,A21t=ω22−μR230ω̇20−μR230−ω̇20ω22+2μR23,A22t=002ω2000−2ω200,B1t=03×3,B2t=I3×3,where μ is the gravitational constant; ω2=a2μ1−e22/R22 is the orbital angular velocity of the target, and a2 is the semimajor axis of the target orbit.For impulse interception, the external acceleration is 0. Letf20 be the true anomaly of the target’s initial position, and the state equation can be solved as follows:
(9)Xtf2=Φtf2,tf20Xtf20,where tf20 is the initial moment, tf2 is the terminal moment, Xtf20 is the relative state of the satellite with respect to the target at the initial moment, Φtf2,tf20 is the state transition matrix, and Xtf2 is the relative state of the satellite with respect to the target at the terminal moment. The division is performed in and out of the orbital plane, and the expression of Φtf2,tf20 is as follows [6]:
(10)Φtf2,tf20=Φ11tf2,tf20Φ12tf2,tf20Φ21tf2,tf20Φ22tf2,tf20=O11tf2,tf200O12tf2,tf2000g1tf2,tf200g2tf2,tf20O21tf2,tf200O22tf2,tf2000g3tf2,tf200g4tf2,tf20.The state transition equation in the orbital plane is(11)where
(12)where k=μ/h3/2, h is the orbital angular momentum of the target, ρf2=1+e2cosf2, ρf20=1+e2cosf20, sf2=ρf2sinf2, sf20=ρf20sinf20, cf2=ρf2cosf2, cf20=ρf20cosf20, sf2′=cosf2+e2cos2f2, cf2′=−sinf2−e2sin2f2, and J=k2tf2−tf20=∫f20f21/ρ2τdτ.The state transfer equation outside the orbital plane is(13)where
(14)g1tf2,tf20g2tf2,tf20g3tf2,tf20g4tf2,tf20=ρf20−e2sinf21k2ρf2−1cosf2−f20sinf2−f20−sinf2−f20cosf2−f20ρf200−e2sinf201k2ρf20.It can be seen from equations (11) and (13) that the state transition matrix of the relative motion between the satellite and the target is converted into a form related only to the initial true anomaly f20 of the target and the terminal true anomaly f2.


## 3.2. Tangent Impulse Interception Solution under the Elliptical Orbit

(15)T=−sinucosΩ−cosucosisinΩ−sinusinΩ+cosucosicosΩcosusini−sinisinΩsinicosΩ−cosi−cosucosΩ+sinucosisinΩ−cosusinΩ−sinucosicosΩ−sinusini.First, the relative position and relative velocity of the satellite with respect to the target in the VVLH coordinate system are calculated. The coordinate transformation matrix from the J2000 coordinate system to VVLH is shown in formula (15), where Ω is the right ascension of ascending node, u is the argument of latitude, and i is the orbit inclination.The coordinate system conversion matrix of the satellite and the target calculated by equation (15) is T1 and T2, respectively. Since the target VVLH coordinate system is a rotating coordinate system, the relative state of the satellite with respect to the target can be expressed in this coordinate system as
(16)r=T2R1−R2,(17)v=T2V1−V2−ω2×r,where R1,V1 and R2,V2 are the positions and velocity vectors of the satellite and the target in the J2000 coordinate system, respectively, and ω2× is the skew-symmetric matrix represented by the angular velocity vector of the target VVLH coordinate system relative to the J2000 coordinate system.
(18)ω2×=00−ω2000ω200.Equation (17) can be further expressed as
(19)v=T2V1−V2+ω2z0−ω2x.Therefore, the relative positionrf20 and the relative velocity vf20 at the initial moment tf20 can be obtained according to equations (16) and (19).Without drifting, the satellite shoots at the initial moment, according to equations (9) and (10); the relative state of the projectile with respect to the target at the terminal moment is
(20)rd,f2vd,f2=Φ11tf2,tf20Φ12tf2,tf20Φ21tf2,tf20Φ22tf2,tf20rd,f20vd,f20,where rd,f2 and vd,f2 are the relative position and velocity vectors of the projectile with respect to the target at the terminal moment, respectively, and the constraint of rd,f2=0,0,0T is satisfied for the interception problem; rd,f20 and vd,f20 are the relative position and velocity vectors of the projectile with respect to the target at the initial shooting moment, and rd,f20=rf20 is satisfied. Therefore, the velocity vector of the projectile relative to the target at the initial moment can be obtained as
(21)vd,f20=−Φ12−1Φ11rf20.Therefore, in the target VVLH coordinate system, the impulse vector launched by the satellite is(22)Δv=vd,f20−vf20.According to formula (19), the velocity vector of the projectile at the moment of shooting in the J2000 coordinate system is
(23)Vd,f20=T2,f20−1vd,f20−ω2,f20zd,f200−ω2,f20xd,f20+V2,f20.Then, the satellite’s launch impulse vector can be expressed in the satellite VVLH coordinate system as(24)Δv1=T1,f20Vd,f20−V1,f20=T1,f20T2,f20−1Δv.And its unit direction vector can be expressed as(25)Δv^1=Δv1Δv12=Δv^1xΔv^1yΔv^1z.Since the orbits of the satellite and the target are coplanar,Δv^1y=0 is satisfied.The flight-direction angle of the initial moment of the satellite is expressed as a univariate function form with respect to the true anomaly of the satellite:(26)γ¯f20=π2−arctane1sinf101+e1cosf10,where f10 is the true anomaly of the satellite at the initial moment.In order to get the true anomaly of the target corresponding to the interception point, define the functionF:
(27)F=−Δv^1x−Δv^1ztanγ¯f20.In the formula, when the initial moment is determined, the flight direction angleγ¯f20 is constant. At this time, the function F is a univariate function about the true anomaly f2 of the target.When the satellite shoots in a tangent mode, the shooting impulse direction of the satellite is consistent with the velocity direction of the satellite at the moment, that is,(28)Ff2=0.Then, according to the change of the true anomaly of the target, the total interception time of the projectile from the initial launch moment to the encounter moment is obtained:(29)ttotal=tf2−tf20=1k2∫f20f21ρ2τdτ.For equation (29), the Lobatto quadrature method can be used for the numerical integration calculation [23].


## 3.3. Solving the Tangent Impulse Intercept Solution

Since using equation (28) is difficult in obtaining an analytical x of the equation is obtained by the combination of incremental search and the false position method [23]. The specific steps are as follows:Step 1.
According to the state of the initial impulse of the satellite, the feasible domain of the interception solution obtained from the second section isf2∈f2,s1,f2,s4.Step 2.
Using the incremental search method, divide the intervalf2,s1,f2,s4 into n adjacent cells, f2,s1,f2,s1,1, f2,s1,1,f2,s1,2,…, f2,s1,n‐1,f2,s4, and determine whether the function symbol of each pair of interval endpoints changes.Step 3.
Use the false position method to obtain an accurate interception solution in the interval where the symbol changes.IfFf2,s1,mFf2,s1,m+1<0, let
(30)f2,r=f2,s1,m+1−Ff2,s1,m+1f2,s1,m−f2,s1,m+1Ff2,s1,m−Ff2,s1,m+1.Takingf2,r as the new endpoint of the interval, it is iterated by equation (30). When f2,r−f2,r+1<δ, the iteration ends, and the exact numerical solution of the interception solution in interval f2,s1,m,f2,s1,m+1 is obtained.


## 4. Minimum Time Tangent Interception Solution Based on GA-SQP

The interception time is different for different impulse points; therefore, the next step is to study the minimum time interception problem in the target orbital period.Let the satellite’s impulse moment bef2t, then the initial drift range is f20→f2t. The relative position and velocity of the satellite with respect to the target at the impulse point can be obtained according to equation (9), namely,
(31)rd,f2tvd,f2t=Φ11tf2t,tf20Φ12tf2t,tf20Φ21tf2t,tf20Φ22tf2t,tf20rd,f20vd,f20.Next, according to Section3, the interception solution corresponding to the satellite at the impulse point f2t can be obtained. For solving the minimum time interception solution (f2∈f2t,f20+2π) in a target period, it takes too much computation to directly calculate the interception time corresponding to each interception solution and find the minimum time solution. Therefore, the paper takes the true anomaly of the target at the satellite’s shooting point and interception point as variables and uses the satellite’s shooting impulse direction, energy, moment, and feasible domain of the solution as constraints to optimize the minimum time interception solution. The formulation of the optimization problem is shown as follows:
(32)minttotal=min1k2∫f2tf21ρ2τdτs.t.Ff2=0Δv12≤maxΔvlimf2t<f2f2t∈f20,f20+2πf2∈f2t,f20+2πf2∈f2,s1,f2,s4.Considering that GA’s global search ability is strong but local search efficiency is low, SQP local search ability is strong but sensitive to initial values. Therefore, the paper combines the advantages of both when solving equation (32). The optimization result of GA is used as the initial value of the SQP search, and the minimum time interception solution in a target period is obtained. The entire algorithm flow is as follows:Step 1.
In population initialization, in the range off20,f20+2π, the satellite’s impulse moment and the interception moment of the projectile are randomly generated and used as the initial population.Step 2.
Determine whether the initial population meets the constraints. Iff2t and f2 generated in Step 1 satisfy the constraints in equation (32), then the fitness function value of the population is calculated. Otherwise, modify the scheme and continue with Step 2.Step 3.
Determine the optimal scheme according to the fitness function value obtained in Step2. If the scheme is optimal and meets the genetic termination condition, the scheme is the output. Otherwise, the program continues to be genetically manipulated.Step 4.
In genetic manipulation, the selection, crossover, and mutation operations are carried out for each group off2t and f2 codes, and the new scheme is obtained by recombining the code segment. The evolutionary generation is set as Q→Q+1. Here, every genetic operation can be considered that the population has evolved a generation. So evolutionary algebra can represent the number of times a population evolves.Step 5.
IfQ≥Qf, output the current optimal scheme and end the operation. Otherwise, continue to execute Step 2. Here, Qf is the total genetic algebra set by the program. Generally, Qf is selected according to the specific simulation task. And the selection needs to ensure that the final output is relatively stable.Step 6.
The optimal solution obtained by the GA algorithm is used as the initial valueg0 of the SQP search, and the QP subproblem is constructed according to equation (32).Step 7.
Solve the QP subproblem, and get the search directiondk and the Lagrange multipliers λk and μk.Step 8.
Ifdk=0, the current result is output as the optimal solution. If dk≠0, the linear search is performed by the golden section method to determine the step size αk and let gk+1=gk+αkdk.Step 9.
If the current result converges, then it is the optimal solution and output it. Otherwise, the Broyden family Davidon-Fletcher-Powell (DFP) revision formula is used to correct the positive definite approximation matrixHk of the Hessen matrix of the Lagrange function [24]. And return to Step 7 for iteration.


## 5. Accurate Interception Solution under High-Precision Extrapolation Model

The method for obtaining the tangent impulse interception solution by the linear model is given in detail in Section3.3. However, the mechanical environment in the actual space is more complicated; if the above model is directly used to get the interception solution, the obtained solution will have a significant error. At the same time, due to the comprehensive influence of the perturbation of gravitational field, nonspherical perturbation, atmospheric drag, solar radiation pressure, ocean tide, and lunar gravitation, it is difficult to obtain an accurate interception solution by analytical means. Therefore, in the following part, the paper iteratively modifies the tangent impulse interception solution obtained under the linear model by combining the optimization method with the high-precision extrapolation model, so as to obtain the high-precision solution of tangent impulse interception.The objective function is defined as the distanceFrttotal,Δv12 from the projectile to the target in the J2000 coordinate system at the interception moment. In the optimization process, the tangent impulse mode is still used. The property of tangent impulse is not changed, but the velocity magnitude of tangent impulse is changed.
(33)Fr=Rd,f2−R2,f22,where ttotal is the interception time of the projectile from the moment of launch to the moment of encounter, Δv12 is the launch speed of the projectile, and Rd,f2 and R2,f2 are the positions of the projectile and the target at the interception moment in the J2000 coordinate system, respectively, both of which are derived from the orbit extrapolation under the high-precision model.It can be seen thatFrttotal,Δv12≥0. Only when the interception time and the impulse velocity take the appropriate value Fr will get a minimum value of 0. Therefore, the problem of solving high-precision tangent impulse interception solution is transformed into an optimization problem.Due to the lack of corresponding gradient information, the Broyden Fletcher Goldfarb Shanno (BFGS) quasi-Newton method is not applicable here, and the intelligent optimization algorithm is relatively time-consuming. Therefore, the Nelder-Mead simplex method is used to solve the problem.The main idea of the Nelder-Mead simplex algorithm is to construct a simple graph withn+1 vertices in an n-dimensional space. Compare the function values on n+1 vertices to determine the search direction of the next step, and operations such as expansion, contraction, shrinkage, and reflection are taken for the simplex. Meanwhile, a new simplex is constructed by replacing the original worst vertex with a better one. Repeat iteration until the optimal solution of the objective function is approximated. The advantage is that it is not sensitive to the initial value and does not need the gradient information. The objective function is not calculated more than twice at each iteration. Therefore, the calculation speed is faster.Since the number of variables in the optimization process is two, the simplex is a geometric object determined by three-point sets. It can be represented by three two-dimensional vectorsp0, p1, and p2, and the determinant
(34)detp0p1p2111≠0is satisfied, wherein the different points pi=ttotal,i,Δv12iT represent the set of variables corresponding to different iteration times. The entire process of the algorithm is as follows:Step 1.
Initialize. During initialization, the vertices constituting the simplex should be linearly independent, that is,(35)detp0p1p2111≠0.
p0 is initialized according to the tangent impulse interception solution obtained by the linear model, and p1 and p2 are arbitrarily selected.Step 2.
Compare the objective function valuesFrp corresponding to the simplex points p0, p1, and p2, and sort out the best point pl, the intermediate point pm, and the worst point ph.Step 3.
Calculate the center pointpc=pl+pm/2 of the simplex points corresponding to two smaller objective functions.Step 4.
Reflect the simplex pointph corresponding to the maximum objective function about pc in the direction of phpc; get the reflection point pr.Step 5.
Calculate the objective function valueFrpr at the reflection point pr, and compare it with Frpl and Frph. According to the possible cases, operations such as expansion, contraction, shrinkage, and reflection are performed to update the simplex [25]. Repeat the iteration until the objective function converges.Based on Sections2–5, summarize the complete process of solving the minimum time tangent impulse interception solution in this paper. The whole procedure is briefly summarized as follows:
(1)
All feasible solutions are obtained within the feasible domain by a numerical iteration method without considering the initial drift segment
(2)
Considering the initial drift segment, the GA-SQP combination optimization algorithm is adopted to obtain the minimum time tangent interception solution within a target period
(3)
The Nelder-Mead simplex method is used to optimize the impulse velocity increment and interception time so as to obtain the accurate interception solution under the high-precision extrapolation model


## 6. Simulation Examples

It is assumed that the flight orbit of the satellite and the target is coplanar, and the length of the earth radiusRe is 6378.13 km. The parameters of the initial flight orbit of the satellite and the target are shown in Table 1.Table 1
Initial orbit parameters of the satellite and target.
Orbital parameter
Target orbit
Satellite orbit
Semimajor axisa (km)
16545.13
16535.13
Eccentricitye
0.4
0.399
Orbit inclinationi (°)
90
90
RAANΩ (°)
0
0
Argument of perigeeϖ (°)
20
19.99
Initial true anomalyf (°)
30
29.99According to the orbital element given in Table1, the relative position and relative velocity of the satellite with respect to the target in the J2000 coordinate system are 9.1596,0,5.3051T km and 4.220,0,4.417T m/s, respectively, at the initial moment. According to equations (16) and (19), the relative relationship between the satellite and the target at the initial moment in the target VVLH coordinate system is rf20=−3.6066,0,−9.9516T km, vf20=−13.024,0,3.1900T m/s, rf202=10.585 km, and vf202=13.409 m/s.According to the initial conditions and the discussion in Section2, the existence interval of the true anomaly of the target at the interception moment is f2,s1,f2,s4=0,2π. It can be seen that the interception solution exists in the target orbit period, and the curve of F in the period is drawn in Figure 2.Figure 2
F value corresponding to the true anomaly of the target under different intercept points.According to Section3.3, the true anomaly of the interception point satisfying condition F=0 is f2=198.6098°, and the impulse energy consumed is Δv12=1.3854 m/s. The unit impulse velocity vector in the J2000 coordinate system is −0.66326,0,0.74839T, which is the same as the flight direction of the satellite, so the shooting was a tangent impulse. The total interception time calculated by equation (29) is 12166.4 s.Since the model used does not take into account the perturbation effect, the situation in a high fidelity model is analyzed next. A HPOP high-precision orbital extrapolation model is selected, which adopts the joint gravity model (JGM4), Jacchia-Roberts atmospheric resistance model, lunar gravity body model, and solar gravity ball model, and takes into account the influence of perturbation of solar radiation pressure, moon tide, and earth tide. The Runge Kutta Fehlberg (RKF7 (8)) integrator is used for the orbital integration, and the integral step is 1 s (the interception time is rounded). In the J2000 coordinate system, the relative position relationship between the projectile and the target obtained under the high-precision model can be seen in Figure3:Figure 3
Under the high-precision model, the distance from the projectile to the target during the interception process in the o-xz plane.The “J” in the coordinate axis in the figure indicates the J2000 coordinate system. It can be seen from the figure that the relative position of the projectile with respect to the target at the interception point is296.449,−0.0585,−67.471T m, and the terminal error is 304.03 m. This error is mainly caused by the model not considering perturbation.In order to obtain the accurate interception solution under the high-precision model, the impulse velocity and interception time obtained under the linear model are taken as the initial values, and the simplex iteration is performed with equation (33) as the objective function. Set the integration step size to 1 s and the termination tolerance to 10−4. The whole iteration process is shown in Figure 4.Figure 4
Convergence process of the relative distance at the terminal moment.Figure4 shows that the terminal distance between the projectile and the target is greatly reduced after the simplex iteration, and the objective function converges after 75 iterations. The relative distance from the projectile to the target reaches −0.0248,−0.0578,−0.0894T m, the terminal error is 0.1093 m, and the impulse energy and the total interception time under the high-precision extrapolation model are Δv12=1.3920 m/s and ttotal=12173 s, respectively.It can be seen that the impulse velocity and interception time obtained under the linear model have high accuracy, which can be used as the initial point for the search of the exact interception solution under the high-precision extrapolation model. After iteration, in the J2000 coordinate system, the relative position relationship between the projectile and the target obtained under the high-precision model can be seen in Figure5.Figure 5
The distance from the projectile to the target during the interception process under the high-precision model after iteration.
(a)
The distance from the projectile to the target during the interception process in theX-axis after iteration
(b)
The distance from the projectile to the target during the interception process in theZ-axis after iterationThe relative positions of the projectile with respect to the target in theY-axis before and after the iteration are the same. This is because when the simplex iteration is performed, the tangent impulse used is located in the o-xz plane, and the velocity compensation of the projectile is not carried out in the Y-axis direction.When the satellite does not shoot at the initial moment, equation (31) can be used to obtain the relative position and velocity of the satellite with respect to the target after drifting. According to Sections 3.3 and 5, the tangent impulse interception solution in period f2t,f2t+2π can be obtained, as shown in Table 2.Table 2
Tangent impulse interception solutions at different impulse points.
f2t (°)
f2 (°)
Two body model
HPOP model
After the simplex iteration
HPOP model
Integration time
Δv12 (m/s)
ttotal (s)
Fr (m)
Δv12 (m/s)
ttotal (s)
Fr (m)
Step (s)
30
198.6098
1.3854
12166
304.03
1.3920
12173
0.1093
1
60
176.7979
1.8140
8643
48.03
1.8131
8662
0.3871
1
90
150.4100
3.4706
4397
153.12
3.4633
4402
0.01948
1
120
471.5169
1.2750
20616
86.30
1.2736
20613
0.8158
1
150
478.1595
1.5628
18405
57.82
1.5381
18392
0.8648
1
180
492.6815
1.8615
15957
370.49
1.8516
15945
0.6607
1
210
508.8669
2.1099
13931
497.95
2.0788
13921
1.0496
1
240
524.0135
2.3256
12992
4745.4
2.1811
12916
1.2098
1
270
272.2454
1361.68
100.13
2.9365
1365.03
100.02
0.2606
0.01
541.7824
2.6133
13486
14719.6
2.2657
13267
1.3267
1
614.9184
3.1056
20417
25974.1
2.7820
20347
1.9056
1
300
304.5269
1085.42
139.26
2.5864
1085.18
139.29
0.3927
0.01
330
335.5802
1126.38
137.55
9.4349
1125.83
137.61
0.3562
0.01
360
366.1213
1051.60
141.61
7.8442
1050.94
141.71
0.5258
0.01In Table2, the satellite transmits a reverse impulse when the shooting moment of the satellite is 270° (f2=272.2454°), 300°, 330°, and 360°. The interception time of the reverse impulse is shorter, the impulse speed is higher, and the integral step has a greater influence on the terminal error. Therefore, the integral step size is shortened to be 0.01 s, and for other shooting moments, it is 1 s. It can also be found that the terminal error under the high-precision model is greatly reduced after the simplex iteration, and the maximum error is 1.9056 m.If the true anomaly of the target is270° when the satellite transmits the impulse, according to the existence condition of the interception solution given in the second section, the feasible range of the true anomaly of the target at the interception point is 0,2π. In the interval 270°,630°, the curve of F is shown in Figure 6.Figure 6
TheF value corresponding to the true anomaly of the target at different interception points.Figure6 shows that if the satellite shoots at this moment, there are three tangent impulse interception solutions. For the impulse points f2t∈f20,f20+2π in a target orbital period, it is required to find out the minimum time interception solution that satisfies the impulse energy requirement.Set the initial population size to 30, the crossover probability to 0.8, the mutation probability to 0.2, the maximum genetic algebra to 20, and the impulse energy to no more than 1.5 km/s. According to equation (32), the minimum time tangent interception solution in a target orbit period is calculated by the GA-SQP algorithm. The iterative process is shown in Figure 7.Figure 7
Iterative process for solving the minimum time tangent interception solution.
(a)
GA iterative process
(b)
SQP iterative processIt can be seen from Figure7 that the GA algorithm converges after the 9th iteration and the minimum time interception solution obtained in the global range is 104.0968 s. The impulse point and the interception point are f2t=271.1563° and f2=273.5302°, respectively. The impulse energy is Δv12=1.3175 m/s. The result obtained by GA is taken as the initial search point of SQP, and the real global minimum time tangent impulse interception solution is obtained after the 4th iteration. The impulse point and the interception point are f2t=267.3343° and f2=269.2662°, respectively. The total interception time is ttotal=89.6197 s. The impulse energy is Δv12=1.5 km/s.In order to verify the accuracy of the results obtained by GA-SQP, the target’s true anomaly at the interception point, velocity increment, and the interception time that are corresponding to each impulse point within the target period is directly calculated at an interval of1°. When there are multiple tangent impulse solutions for one impulse point, the time minimum solution is selected. In this way, for a fixed impulse point f2t∈f20,f20+2π in the target orbit period, there is only one corresponding interception point. Therefore, the interception solutions of tangent impulse at different impulse points in the target period can be obtained. The result is shown in Figure 8.Figure 8
Tangent impulse interception solution at different impulse points in the target orbit period.
(a)
The target’s true anomaly at the interception point corresponding to different impulse points
(b)
Velocity increment at different impulse points
(c)
Interception time corresponding to different impulse pointsAccording to Figures8(a) and 8(b), the impulse points satisfying the constraint conditions are f2t∈267°,383°∪30°,112° within the target orbit period. Meanwhile, as shown in Figure 8(c), the impulse point corresponding to the minimum time tangent impulse interception solution is f2t∈267°,268°. In summary, it can be seen that the solution obtained by the GA-SQP algorithm is very close to the global minimum time tangent impulse interception solution.Taking the integral step to 0.01 s, the simplex method is used to iterate out the minimum time tangent impulse interception solution under the HPOP model. After the iteration, the relative position of the projectile with respect to the target in the J2000 coordinate system is−0.1738,5.6011×10−5,−0.4076T m, and the terminal error is 0.4431 m. Moreover, the impulse energy and the interception time are Δv12=1.5 km/s and ttotal=89.62 s, respectively. The above results show that the minimum time tangent impulse interception solution calculated by the linear model has high accuracy, and it does not change substantially before and after the iteration. The entire interception process is shown in Figure 9.Figure 9
Interception process in the J2000 coordinate system.


## 7. Conclusion

The paper studies the minimum time interception problem under tangent impulse by using the relative motion equation of the elliptic orbit. The problem is transformed into a univariate function for solving the true anomaly of the target by a linear relative motion equation. When the impulse point is given, the feasible domain of the target true anomaly is solved firstly. Then all feasible solutions are obtained by incremental search and a trial position method. When the impulse point is freely selected, the minimum time tangent impulse interception solution is obtained by the GA-SQP combination algorithm. Considering the perturbation and nonlinear effects, the interception solution is optimized by the Nelder-Mead simplex method so as to obtain the accurate interception solution under the high-precision orbit extrapolation model.Simulation examples show the following:(1)
This method is suitable for the case that the target orbit is an ellipse. The obtained interception solution has high accuracy and can be used as an initial point for the search of the exact interception solution under the high-precision extrapolation model
(2)
When gradient information is not easily accessible, the simplex iteration can greatly reduce the satellite’s shooting terminal error under the high-precision extrapolation model. However, this conclusion is only applicable to the case that the target is relatively close to the satellite. When the two are relatively far away, the optimization performance of simplex iteration may be degraded due to the decrease of the initial solution accuracy
(3)
The GA-SQP combination algorithm effectively combines the advantages of GA and SQP. It not only overcomes the weak local search ability of GA but also avoids the problem that SQP is sensitive to the initial value. The exact minimum time tangent impulse interception solution can be obtained by using this algorithm


### 7.1. Prospects for Future Work

Due to the influence of nonlinearity and perturbation, when the relative distance between the target and the satellite is far away, the accuracy of the initial solution provided by the model will get worse. Therefore, the next step can be to analyze the effect of the change in relative distance between the target and the satellite on the optimization effect of the simplex.


## 7.1. Prospects for Future Work

Due to the influence of nonlinearity and perturbation, when the relative distance between the target and the satellite is far away, the accuracy of the initial solution provided by the model will get worse. Therefore, the next step can be to analyze the effect of the change in relative distance between the target and the satellite on the optimization effect of the simplex.



---

*Source: 6827982-2019-08-05.xml*