Research Article  Open Access
Optimal Midcourse Guidance Algorithm for Exoatmospheric Interception Using Analytical Gradients
Abstract
This paper is aimed at providing a semianalytical method to solve the optimal exoatmospheric interception problem with the minimum fuel consumption. A nonlinear programming (NLP) problem with the minimum velocity increment, which involves Lambert’s problem with unspecified timeofflight, is firstly formulated. Then, a set of KarushKuhnTucker conditions and the Jacobian matrix corresponding to those conditions are derived in an analytical manner, even though the derivatives are mathematically complicated and computationally onerous. Therefore, the NewtonRaphson method can be used to efficiently solve this problem. To further decrease computational cost, a nearoptimal initialization method reducing the dimension of the search space is presented to provide a better initial guess. The performance of the proposed method is assessed by numerical experiments and comparison with other methods. The results show that this method is not only of high computational efficiency and accuracy but also applicable to onboard guidance.
1. Introduction
Exoatmospheric missile defense system is a generic term conveying an interceptor designed to destroy any ballistic targets delivering nuclear, chemical, biological, or conventional warheads outside the atmosphere. It is of great significance for national security, and therefore, the majority of powerful countries such as the United States, Russia, and China have developed their own missile defense systems [1–3]. In the midcourse, both interceptor and target follow Keplerian orbits and are in highspeed interception engagement with the relative velocity greater than 10 km/s. Thus, “hittokill” is the only valid way of completely destroying invading targets, which brings forward high requirements on the performance of guidance system [4–6]. Actually, exoatmospheric midcourse guidance is a process to calculate an initial velocity vector to ensure a successful impact with zero miss, which theoretically can be categorized into a type of Lambert’ problem [7–9]. Our purpose in this paper is to develop a guidance algorithm from the point of solving Lambert’s problem efficiently.
In 1809, Gauss [10] developed the first iterative process to solve Lambert’s problem, but it suffers from the singularity problem. Battin et al. [11, 12] improved Gauss’s method and proposed an elegant algorithm that efficaciously removed the singularity. Their method is based on a new transformation and a new iteration function to achieve fast convergence. It is mathematically elegant and practically implementable, but its derivation is founded on the complicated geometric properties of conic sections. Avanzini [13] presented an intuitive and simplified version of Battin’s method by parametrizing admissible orbits in terms of the transverse eccentricity component. An iterative method based on Householder’s root solver was designed in [14], where the first and second derivatives of the timeofflight (TOF) equation are analytically derived to increase the rate of convergence. In [15], LeviCivita regularization was applied to Lambert’s problem, and a NewtonRaphson method in combination with safety checks was adopted to achieve both speed and robustness. These works have successfully promoted the application of Lambert’s problem in the space science community. Additionally, Nelson and Zarchan’s efficient method [16, 17] was popular in defense applications [18]. This method regarded the flight path angle as an iterative variable and employed a secant method combined with boundary values to speed up the routine. Obviously, the aforementioned researchers focused on numerical procedures to deal with Lambert’s problem by searching a variable iteratively, because it is impossible to derive an analytical solution due to transcendental equations. Therefore, it is still of great significance to speed up the convergence of iterative algorithms, especially for onboard autonomous scenarios requiring computation [19–22].
Moreover, for some situations where the issue of convergence time is critical, two types of methods are commonly used to improve convergence rate. The first one focuses on analytical gradients [23]. Bate et al. [24] presented an effective algorithm named as the iteration method, where the slope of the TOF with respect to the semilatus rectum could be expressed in an analytical manner. Ahn and Lee [25] analytically derived the gradient of the flight path angle with respect to the TOF from the conditions of the twopoint orbital boundary value problem, and their method was fully verified by comprehensive numerical experiments. Furthermore, another practical way to accelerate the convergence of iteration process is improving the initial guess. In 1990, Gooding [26] made pioneering efforts in this field and used Halley’s iterations to set the initial guess, thereby achieving fast convergence properties. Then, Izzo [27] improved Gooding’s method by inverting the linear approximation to the TOF curves, which significantly reduced the computational complexity. Arora and Russell [28] proposed a fast and robust algorithm to solve multiple revolution Lambert’s problem, in which a new geometry parameter was introduced to simplify the TOF equation, and therefore, an accurate initial guess can be provided for rapid rootsolving. In [29], Ahn and Lee also developed a method using a twodimensional table involving the geometric characteristic and the normalized TOF to interpolate an initial guess, whose efficiency was validated by various numerical experiments. It should be noted that the aforementioned methods assume that the TOF is required to be a given constant.
The optimal exoatmospheric interception problem with the minimum fuel consumption is to find the optimal velocity increment to eliminate interception error and retain enough energy for a successful collision. This kind of problem is abbreviated as the minimum velocity increment problem (MVIP) in this paper. It is worth noting that the aforementioned algorithms cannot be directly applied to the MVIP, because the TOF should be regarded as an unspecified parameter in the optimization process [30–32]. Commonly, an additional singlevariable algorithm needs to be developed to repeatedly solve Lambert’s problem so as to determine the optimal TOF. However, this type of method has high computational cost because of its overreliance on a twolevel iterative process. For onboard autonomous scenarios, the computation time for generating commands is very critical. Therefore, it is necessary to develop a new algorithm for solving the MVIP in an inexpensive manner.
The core objective of this paper is to provide a semianalytical method for solving the MVIP using analytical gradients, and it has potential advantage in onboard guidance. The main contributions of this paper are stated as follows. Firstly, a typically nonlinear programming problem (NLP) with the minimum velocity increment is formulated, in which two transcendental equations characterizing Lambert’s problem and Keplerian motion are regarded as equality constraints in order to reduce computational difficulties. Secondly, a set of KarushKuhnTucker (KKT) conditions for this NLP problem is derived to determine the minimum velocity increment. To speed up convergence, the Jacobian matrix is successfully derived on each parameter in an analytical way, even though the expressions are mathematically complicated and computationally onerous. Accordingly, rootfinding algorithms such as the NewtonRaphson method can be effectively used to find the optimal solution with a high accuracy. Thirdly, it is found that the variable, semilatus rectum , is insensitive to the variation of TOF. A nearoptimal initialization method, which considers as a constant and reduces the dimension of the search space, is developed to provide a better initial guess. That further accelerates the rate of convergence. In a word, the proposed method is attractive from the point of view of low computational cost, high accuracy, and the fact that it has potential to be used as a baseline algorithm for onboard guidance. Finally, the performance of this method is sufficiently verified by comprehensive experiments and comparison with the existing methods.
The rest of the paper is organized as follows. Section 2 provides a detailed description of the MVIP, which involves orbital motion and Lambert’s problem, respectively. Section 3 introduces the method used to solve the problem in this paper, which is composed of the KKT conditions, analytical gradients, initialization method, and the procedure of implementation. Section 4 presents the numerical results to demonstrate the performance of the method. Finally, conclusions are given in Section 5.
2. Problem Formulation
The engagement geometry of the MVIP is illustrated in Figure 1(a). At the moment , the target performs an orbital maneuver and its trajectory is transferred from the Keplerian orbit orb1 to a new orbit orb2. The exoatmospheric kill vehicle (EKV), which is a representative interceptor designed for exoatmospheric interception, needs to change its velocity vector from to so as to score a direct hit on the predicted intercept point (PIP) PIP_{2}. Commonly, a PIP can be uniquely determined by orbital equations if the TOF is given. Then, the midcourse guidance problem can be formulated as the determination of an orbit having a specified transfer time and connecting two position vectors, which is known as Lambert’s problem. Actually, there must be an optimal TOF corresponding to the minimum velocity increment , where the symbol “” denotes the Euclidean norm of vector. To demonstrate the MVIP more clearly, a schematic representation for the correlation between Keplerian motion and Lambert’s problem is presented in Figure 1(b).
(a)
(b)
2.1. Orbital Motion of the Target
Let the orbital elements of the target orbit orb2 be where is the right ascension of the ascending node, is the inclination of the orbit, is the semimajor axis, is the eccentricity, is the argument of periapsis, and denotes the time of periapsis passage. The subscript denotes the invading target.
The true anomaly is defined as the angle from the periapsis to the orbiting object measured in the direction of motion, which is denoted as . The target position vector, , and the corresponding geocentric distance, , are expressed as follows: where denotes the unit vector of terminal position and is defined as follows:
Furthermore, the geometric relationship between and its corresponding eccentric anomaly is and according to orbital motion equations, the transfer time, , can be formulated as follows:
Then, substituting equation (5) into equation (6), the orbital motion of the target can be explicitly expressed as a function of and .
Equations (2) and (7) indicate that the position vector, , is uniquely determined by the TOF. However, equation (7) is transcendental in and unable to be solved analytically.
2.2. Lambert’s Problem
As stated above, if the terminal position vector of the target, , is determined, the midcourse guidance problem can be formulated as a classical Lambert’s problem. Its solution is the initial velocity vector, which ensures the interceptor will impact the target with zero miss. As the Keplerian orbit is confined to a plane, Lagrangian coefficients can be used to express the equations for terminal vectors, and , in terms of initial vectors, and ; thus,
Clearly, the coefficients, and , are simply the derivatives of and with respect to the time, respectively, and they can be expressed in terms of cross products as follows: where is the magnitude of , which is the angular momentum vector of the EKV.
Furthermore, the position and the velocity vectors are expressed in terms of orbital plane coordinates as follows: where is a unit vector in the direction of perigee, the point of closest approach to the earth, and is in the orbital plane and normal to such that complete a righthanded coordinate system The equation of an orbit in polar coordinates is written as follows: where is the semilatus rectum, and denote the semimajor axis, the eccentricity, and the true anomaly of the EKV orbit, respectively. Differentiating equation (12) and one can have where is the standard gravitational parameter, which equals 3.986 × 10^{14} m^{3}/s^{2}. Then, substituting equations (10) and (13) into equation (9), the Lagrangian coefficients can be rewritten as follows: where denotes the angle between and ; thus,
It should be noted that, s a constant in the classical Lambert’s problem, and therefore, the Lagrangian coefficients are uniquely determined by .
In addition, the position and velocity vectors in equation (10) can be directly written in terms of the eccentric anomaly as follows:
Taking the derivative of equation (6), one can obtain that
Substituting equations (17) and (19) into equation (9), and the Lagrangian coefficients can be rewritten as follows: where denotes the difference of the eccentric anomaly, which is expressed as .
Combining equations (15) and (20), the transfer time can be expressed as a function of Lagrangian coefficients where the semimajor axis is expressed as follows:
Substituting equations (22) and (15) into equation (21), Lambert’s problem can be finally expressed as a function of and
It is apparent that equation (23) is a transcendental equation, and therefore, its solution can only be provided by iterative methods.
In summary, (equation (7)) and (equation (23)) characterize the numerical properties of Keplerian motion and Lambert’s problem, respectively. Once a TOF is determined, the intercept orbit can be ascertained by solving equation (23), in which equation (7) is an equal constraint used to determine the terminal positon vector . Then, substituting the value of into equation (15), the initial velocity vector of the EKV can be calculated by Lagrangian coefficients as follows:
Because is uniquely defined by the TOF, a singlevariable unconstrained NLP problem can be formulated for the purpose of determining the minimum velocity increment; thus, the cost function is defined as follows: where the square of the velocity increment is given by
To find a local minimum , it is necessary to derive a set of firstorder necessary conditions of this problem
It should be noted that two transcendental equations, (equation (7)) and (equation (23)), are involved in the computational process, which leads to the strong coupling among if the transfer time is considered as an only independent variable. Therefore, the mathematical complexity is sharply increased. In order to overcome the computational difficulties yielding from transcendental equations, the independent variable is expanded to . Meanwhile, are considered as two equality constraints. Although the number of iterative variables increases, additional items resulting from the coupling are efficiently eliminated. Subsequently, the NLP problem formulated in equation (25) is transferred to a multivariable NLP program with equality constraints. The exhaustive derivations of KKT conditions and analytical gradients are presented in Section 3.
3. KKT Condition Solution Procedure Using Analytical Gradients
This section introduces the whole procedure to solve the MVIP based on the numerical computation of the NLP problem. Section 3.1 presents the mathematical basis of NLP and obtains the KKT conditions corresponding to the MVIP. Then, the NewtonRaphson method is introduced to obtain the optimal guidance command. Section 3.2 derives the Jacobian matrix analytically on each parameter, which increases the rate of convergence of the NewtonRaphson iteration. Section 3.3 develops an initialization method used to conduct a nearoptimal initial guess for further accelerating the convergence. Section 3.4 describes the strategy to update the normalized iterative variables.
3.1. Derivation of the KKT Conditions
Consider a NLP problem with both equality and inequality constraints. The objective is to minimize the cost function subject to a set of equality constraints and inequality constraints . The problem is then
The KKT conditions are found by defining the augmented cost function, , so that where are the KKT multipliers associated with equality and inequality constraints, respectively. The optimal solution, , is determined by the gradient of the augmented cost function [33]. If are continuously differentiable functions, and then there exists KKT multipliers and such that
For the MVIP discussed in this paper, the corresponding NLP problem is finally expressed as follows:
It should be emphasized that inequality constraints can be omitted when the iterative variables are restricted to the definition domain. In addition, the terminal position of Lambert’s problem is no longer a given constant in the MVIP, but is determined on the orbit of the target as formulated in . Thus, is involved in the third equation of equation (33), which also decreases the computational difficulty. Therefore, the independent variables in the equality constraint are set to be , which is different from that in equation (23). Furthermore, the augmented cost function is given by
Taking the partial derivative of the independent variables, one can obtain
Rearranging equation (35), then the KKT multipliers and can be formulated as follows:
Substituting equation (36) into equation (35) and comparing with equation (30), the KKT conditions can be denoted by a system of transcendental equations as follows: where where and describe the mathematical models of Keplerian motion and Lambert’s problem in the MVIP. Meanwhile, relates the minimum velocity increment of the EKV with the orbital motion of the target through the KKT conditions. To solve the system of transcendental equations in equation (38), a NewtonRaphson iterative method is designed as follows.
The improvements of iterative variables are formulated as follows: where the superscript means the th iteration, and the Jacobian matrix is given by
In spite of mathematical complication, analytical expressions of all the parameters in the are precisely derived in the following subsection, thereby significantly accelerating the computational efficiency and improving the accuracy in comparison with finite difference.
3.2. Derivation of the Analytical Gradients
This subsection focuses on the analytical expression of the Jacobian matrix on each parameter. The partial derivatives of are firstly formulated due to its concise structure.
3.2.1. Partial Derivatives of
Since is not explicit in , the partial derivative of with respect to is zero.
The partial derivative of with respect to can be obviously expressed as follows:
Substituting equation (5) into equation (39), the partial derivative of with respect to is given according to trigonometry formulas
Although the partial derivatives of can be formulated succinctly, other formulas discussed later are very complex. Nevertheless, an exhaustive derivation of analytical gradients will be presented.
3.2.2. Partial Derivatives of
Observing the expression of in equation (39), the partial derivative with respect to t can be obviously expressed as follows:
For simplifying the derivations, three notations are introduced as follows:
Substituting equation (15) into equation (22), the semimajor axis of the intercept orbit is rewritten as follows:
Then, the partial derivative of with respect to is expressed as follows: where is another notation expressed as follows:
Substituting equations (15) and (49) into equation (39) and partially differentiating with respect to yields where
Partially differentiating equation (2) with respect to , one can obtain that
Substituting equation (16) into the expression of in equation (15), it can be rewritten as follows:
Partially differentiating equation (55) with respect to, and introducing notations as follows: where
Then, the partial derivative of the Lagrangian coefficient with respect to is given by
In order to simplify the expression, notations and are also introduced, which are expressed as follows:
Partially differentiating equation (49) with respect to , and substituting equation (54) into the result, one finds
Then, substituting equations (15) and (55) into equation (39) and partially differentiating with respect to yields where
So far, the partial derivatives of have been successfully expressed by equations (47), (52), and (61) in an analytical manner. It is worth noting that consists of the partial derivatives of . Therefore, the partial derivatives of involve a number of secondorder partial derivatives, which are much more complicated.
3.2.3. Partial Derivatives of
Observing the expression for in equation (39), it is easy to see that does not appear explicitly, so the partial derivative of with respect to can be expressed as follows:
Then, let us focus on the partial derivative with respect to . The derivation is presented as follows:
At the beginning, substituting equations (2) and (24) into equation (33) and introducing a set of notations then the cost function can be rewritten as follows:
Partially differentiating equation (65) with respect to and introducing another set of notations as one has
Then, differentiating equations (50), (67), and (52) with respect to, one can obtain secondorder partial derivatives where
Partially differentiating in equation (15) and using equations (16) and (54), one has
Then, differentiating equations (50), (65), and (52) with respect to and introducing two sets of notations expressed as one can obtain the secondorder partial and mixed derivatives
Substituting equations (46), (52), (61), and (67) into equation (39), and partially differentiating with respect to , the partial derivative of with respect to is finally derived as follows:
Let us focus on the last term of the partial derivative of . Partially differentiating equations (54), (58) ,and (70) with respect to , one can obtain that
Partially differentiating equation (65) with respect to twice and introducing another set of notations are expressed as follows: where
Then, the secondorder partial derivative of the cost function with respect to is given by
Partially differentiating equation (60) with respect to then one has where
Partially differentiating equation (61), one can obtain the secondorder partial derivative of with respect to where
Furthermore, partially differentiating equation (46), the secondorder partial derivative of is expressed as follows:
Then, the partial derivative of g_{3} with respect to f_{M} can be expressed as follows:
Now, the partial derivatives of have been successfully presented by equations (63), (73), and (83); thus, all the parameters of the Jacobian matrix have been analytically formulated. Although the derivations are mathematically complicated and computationally onerous, it only takes small fractions of a second to compute the gradient information, which is of great significance to the improvement of the convergence rate. Moreover, to further improve the speed of convergence, a nearoptimal initialization method is developed in the next subsection.
3.3. Initialization Method
In order to accelerate convergence, this section presents an efficient method to provide an improved initial guess of iterative variables for the MVIP. This method is based on a discovery that semilatus rectum, , is insensitive to the variation of TOF in the computation process. In fact, the semilatus rectum is a function of the semimajor axis and the eccentricity. It is impossible for EKV to substantially change these two orbital elements due to the limitation of fuel. Therefore, the semilatus rectum varies within a limited range. Furthermore, an initialization method considering the semilatus rectum as a constant can be used to provide a warm start. The concrete computational procedure is given as follows.
At the beginning, is calculated from initial conditions and considered as a constant in computation process. Therefore, the dimension of the search space is reduced to two, which are the transfer time, , and the true anomaly, . It should be noted that, according to equation (7), can be analytically expressed as a function of . Therefore, the dimension of the search space can be further reduced to formulate a singlevariable rooting problem. Let denote the equality constraint on which is expressed as follows:
To determine the initial guess of a singlevariable NewtonRaphson iteration equation is formulated as follows: where
It should be pointed out that can be calculated analytically from the derivation in Section 3.2.
The iteration is finished if the difference between and is small enough, and the convergence criterion of the initialization method is given by
The proposed initialization method is derived based on the physical characteristic of the exoatmospheric interception. It can not only reduce the difficulty in evaluation by compressing the search space to one parameter but it also provides a feasible and nearoptimal initial guess for the main iteration, thereby efficiently accelerating the rate of convergence.
3.4. Updating Iterative Variables
When the initial guess has been determined, optimal parameters corresponding to the minimum velocity increment can be efficiently provided by the NewtonRapson method with analytical gradients. If the difference between two successive values is smaller than the convergence criterion the main iteration stops and the velocity increment of the current iteration is the optimal solution of the MVIP. The convergence criterion is expressed as follows:
Otherwise, the iterative variables should be updated to start the next iteration.
It should be pointed out that the MVIP has significant magnitude difference among the iterative variables. For example, the magnitude of semilatus rectum reaches as high as 10^{6}, whereas the true anomaly locates on the interval . Therefore, it is essential to normalize the iterative variables to provide a better numerical condition for solving the NLP problem. The normalization is achieved by multiplying a special weight matrix and the corresponding convergence criterion is formulated as follows: where is the weight matrix.
The procedure of implementing the proposed algorithm is included in the flowchart in Figure 2. After the input of program data, the initialization method is applied to provide a warm start for the main iteration. Then, analytical gradients are used to calculate the Jacobian matrix in the computational process of the NewtonRapson iteration, which is propitious for the computational efficiency. It should be noted that the initialization method is only applied at the first step of onboard guidance. The reason is that the optimal result generated from the previous guidance cycle is very close to the optimal one in current cycle and can be considered as the initial guess for reducing computational cost. Furthermore, a comprehensive study on the charming properties of the algorithm is presented in Section 4.
4. Validation of the Proposed Method
This section presents the numerical results used to verify the proposed algorithm. The validation is carried out from three different perspectives. Firstly, the convergence property of the algorithm is investigated for nine cases with different initial conditions. Secondly, comprehensive experiments are conducted to assess the computational efficiency by comparing the number of iteration (NOI) and computation time with other typical methods. Thirdly, the proposed method is implemented in the framework of midcourse guidance, which is used to validate its applicability. The corresponding comparison of fuel consumption with other guidance laws is also provided. All algorithms are implemented in MATLAB and run on personal computer with Intel Corei7 CPU and 16G RAM under the Win10 operating system. It should be noted that all the numerical experiments have different initial conditions which are divided into three groups with different target orbits. Each experiment group corresponds to each of three different initial velocity vectors, which is challenging enough to evaluate the performance of the proposed method. Nine simulation conditions are listed in Table 1 in detail.

4.1. Convergent Characteristics
The three groups of simulation cases listed in Table 1 are used to demonstrate the convergent characteristics of the proposed method including the NOI, the convergent process of iterative variables, and the magnitude of convergence criterion. These three groups have a large difference in their initial heading angles which represent typical trajectories of intercontinental ballistic missiles (ICBMs). Meanwhile, the interceptor in each group has different initial velocity vector, which is designed to ensure that the optimal impact happens at the beginning, middle, and end phase of exoatmospheric ballistic trajectory. The engagement geometry and the iterative process of the cost function are illustrated in Figures 3–5. Obviously, it only takes a few cycles for the proposed method to successfully converge to the optimal results in all experiments. In fact, benefiting from the analytical gradients, the cost function has approached a quite small neighborhood of the optimal solution after the first two iterations. Because of the limitation of the length of paper, Table 2 only presents the summary of the convergence properties for the cases of Group 1.
(a)
(b)
(a)
(b)
(a)
(b)
