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.

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 [1215]. 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 in Figure 1, there is a solution only when the interception point is located at .

In the above figure, the orbits of the target and the satellite are both anticlockwise. is the shooting point of the satellite; v is the velocity of the satellite for this moment; and are the two intersections of the satellite’s velocity direction and the target orbit, respectively; is the intercepting point on the target orbit; is the focus of the target orbit as well as the focus of the satellite orbit; and are radius vectors of the and 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 where where and are the arguments of perigee of the satellite and the target orbit, respectively, and are the true anomalies of the satellite and the target orbit, respectively, and are the eccentricity ratio and semilatus rectum of the target orbit, respectively, and is the distance between the satellite and the center of the earth. (1)If , then (2)If , there is no solution for (3)If , then where

When the transfer orbit of the projectile is an ellipse, the constraint is [9] where where is the transfer angle of the projectile.

As the analysis of equation (1), there are three situations to discuss: (1)If , then (2)If , there is no solution for (3)If , then where

It can be seen from the above conclusion that the feasible domain of the tangent interception solution corresponding to the satellite impulse point is , in which, when , 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 .

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 its -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 -axis points to the direction of motion of the spacecraft, the -axis points to the center of the earth, and the -axis meets the right-hand rule.

In the VVLH coordinate system, the relative motion of the spacecraft with respect to the target is [6] where and ; is the external acceleration, and the coefficient matrix is where is the gravitational constant; is the orbital angular velocity of the target, and is the semimajor axis of the target orbit.

For impulse interception, the external acceleration is 0. Let be the true anomaly of the target’s initial position, and the state equation can be solved as follows: where is the initial moment, is the terminal moment, is the relative state of the satellite with respect to the target at the initial moment, is the state transition matrix, and 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 is as follows [6]:

The state transition equation in the orbital plane is where where , is the orbital angular momentum of the target, , , , , , , , , and .

The state transfer equation outside the orbital plane is where

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 of the target and the terminal true anomaly .

3.2. Tangent Impulse Interception Solution under the Elliptical Orbit

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, is the argument of latitude, and is the orbit inclination.

The coordinate system conversion matrix of the satellite and the target calculated by equation (15) is and , 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 where and are the positions and velocity vectors of the satellite and the target in the J2000 coordinate system, respectively, and is the skew-symmetric matrix represented by the angular velocity vector of the target VVLH coordinate system relative to the J2000 coordinate system.

Equation (17) can be further expressed as

Therefore, the relative position and the relative velocity at the initial moment 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 where and are the relative position and velocity vectors of the projectile with respect to the target at the terminal moment, respectively, and the constraint of is satisfied for the interception problem; and are the relative position and velocity vectors of the projectile with respect to the target at the initial shooting moment, and is satisfied. Therefore, the velocity vector of the projectile relative to the target at the initial moment can be obtained as

Therefore, in the target VVLH coordinate system, the impulse vector launched by the satellite is

According to formula (19), the velocity vector of the projectile at the moment of shooting in the J2000 coordinate system is

Then, the satellite’s launch impulse vector can be expressed in the satellite VVLH coordinate system as

And its unit direction vector can be expressed as

Since the orbits of the satellite and the target are coplanar, 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: where 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 function :

In the formula, when the initial moment is determined, the flight direction angle is constant. At this time, the function is a univariate function about the true anomaly 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,

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:

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 is .

Step 2. Using the incremental search method, divide the interval into adjacent cells, , ,…, , 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.

If , let

Taking as the new endpoint of the interval, it is iterated by equation (30). When , the iteration ends, and the exact numerical solution of the interception solution in interval 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 be , then the initial drift range is . 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,

Next, according to Section 3, the interception solution corresponding to the satellite at the impulse point can be obtained. For solving the minimum time interception solution () 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:

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 of , 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. If and 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 Step 2. 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 of and codes, and the new scheme is obtained by recombining the code segment. The evolutionary generation is set as . 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. If , output the current optimal scheme and end the operation. Otherwise, continue to execute Step 2. Here, is the total genetic algebra set by the program. Generally, 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 value of the SQP search, and the QP subproblem is constructed according to equation (32).

Step 7. Solve the QP subproblem, and get the search direction and the Lagrange multipliers and .

Step 8. If , the current result is output as the optimal solution. If , the linear search is performed by the golden section method to determine the step size and let .

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 matrix 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 Section 3.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 distance 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. where is the interception time of the projectile from the moment of launch to the moment of encounter, is the launch speed of the projectile, and and 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 that . Only when the interception time and the impulse velocity take the appropriate value 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 with vertices in an -dimensional space. Compare the function values on 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 vectors , , and , and the determinant is satisfied, wherein the different points 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, is initialized according to the tangent impulse interception solution obtained by the linear model, and and are arbitrarily selected.

Step 2. Compare the objective function values corresponding to the simplex points , , and , and sort out the best point , the intermediate point , and the worst point .

Step 3. Calculate the center point of the simplex points corresponding to two smaller objective functions.

Step 4. Reflect the simplex point corresponding to the maximum objective function about in the direction of ; get the reflection point .

Step 5. Calculate the objective function value at the reflection point , and compare it with and . 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 Sections 25, 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 radius is 6378.13 km. The parameters of the initial flight orbit of the satellite and the target are shown in Table 1.

According to the orbital element given in Table 1, the relative position and relative velocity of the satellite with respect to the target in the J2000 coordinate system are km and 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