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 time-of-flight, is firstly formulated. Then, a set of Karush-Kuhn-Tucker 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 Newton-Raphson method can be used to efficiently solve this problem. To further decrease computational cost, a near-optimal 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 high-speed interception engagement with the relative velocity greater than 10 km/s. Thus, “hit-to-kill” 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 time-of-flight (TOF) equation are analytically derived to increase the rate of convergence. In [15], Levi-Civita regularization was applied to Lambert’s problem, and a Newton-Raphson 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 two-point 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 root-solving. In [29], Ahn and Lee also developed a method using a two-dimensional 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 single-variable 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 two-level 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 Karush-Kuhn-Tucker (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, root-finding algorithms such as the Newton-Raphson 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 near-optimal 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) PIP2. 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 right-handed 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 × 1014 m3/s2. 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 single-variable 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 first-order 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 Newton-Raphson 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 Newton-Raphson iteration. Section 3.3 develops an initialization method used to conduct a near-optimal 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 Newton-Raphson 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 second-order 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 second-order 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 second-order 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 second-order 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 second-order partial derivative of with respect to where
Furthermore, partially differentiating equation (46), the second-order partial derivative of is expressed as follows:
Then, the partial derivative of g3 with respect to fM 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 near-optimal 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 single-variable rooting problem. Let denote the equality constraint on which is expressed as follows:
To determine the initial guess of a single-variable Newton-Raphson 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 near-optimal 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 Newton-Rapson 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 106, 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 Newton-Rapson 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 Core-i7 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)
As stated above, the proposed method contains two parts: the initialization method, which is used to provide a near-optimal initial guess, and the main iteration which uses analytical gradients to solve the MVIP so as to obtain the minimum velocity increment. As shown in Table 2, two parts successfully converge to the required tolerance within six cycles. The maximal error for equality constraints, , is less than . The error norm of iterative parameters monotonously decreases and its maximal error is . Therefore, the proposed method works well for all cases and has a high convergence rate.
In addition, let us focus on the convergent process of listed in the right column of Table 2, the maximal relative error between the initial and optimal values is barely 14.25% for Case 1, 2.17% for Case 2, and 0.22% for Case 3. Therefore, it is reasonable to assume that can be considered as a constant in the initialization method. Another attention should be paid is that a relaxed convergence criterion is used for the initialization method, which can further reduce the computation time. However, a strict convergence criterion should be employed in the main iteration for obtaining high accuracy results, so that
In summary, the proposed method is able to achieve extremely high accuracy with few iterations and is worthy of practical application. Moreover, the computational efficiency of the proposed method is discussed comprehensively in the next subsection by comparing with other algorithms.
4.2. Computational Efficiency
In order to assess the computational efficiency of the proposed method, a comparison with typical methods is provided in this subsection. The interception scenarios are the same as that listed in Table 1. The methods used for comparison include Zarchan’s method and the -iteration method, which are both widely used in solving Lambert’s problem. The major difference is that Zarchan’s method is based on finite difference, whereas the -iteration method has the analytical gradient. Additionally, the active set method, which is a general algorithm designed to solve NLP problem, is also used for comparison. The numerical results of the NOI and the computation time for different methods are presented in Tables 3 and 4, respectively.
It can be seen from Table 3 that the NOI of the proposed method is only six, which is two to three orders of magnitude smaller than that of Zarchan’s method and the -iteration method, respectively. It should be noted that the two methods are not pertinently designed for the MVIP, and therefore, an additional single-variable algorithm iterating on the TOF should be developed to recurrently solve Lambert’s problem, which leads to high computational burden. However, the active set algorithm is used to solved the NLP problem transferred from the MVIP. The maximal NOI is 30, which is much smaller than that of Zarchan’s method and the P-iteration method. Moreover, the active set algorithm uses finite difference approach to provide gradient information, which leads to worse numerical accuracy than analytical gradients. Therefore, the NOI of the proposed method is two and three times smaller than that of the active set algorithm.
Let us focus on the computation time as shown in Table 4. The average value is barely 0.089 ms for proposed method, but is 7.614 ms for Zarchan’s method, 4.451 ms for the -iteration method, and 19.42 ms for the active set algorithm. It is obvious that the computation time of the proposed method is much smaller than that of the others. It should be pointed out that the active set algorithm costs much time due to the invocation of toolbox functions in MATLAB, whereas the proposed method is independent and only consists of basic operation.
In conclusion, the average values of NOI and computation time are both improved by two to three orders of magnitude, which makes the proposed method to be particularly suited for onboard applications. The corresponding statistic is shown in Figure 6. It should be noted that Zarchan’s method and the -iteration method have been demonstrated to be effective in solving Lambert’s problem, and these reductions do not necessarily mean that the proposed method outperforms that of the conventional algorithms. In fact, the computational challenge of the MVIP arises from the combination of Lambert’s problem and orbital motion, which brings more transcendental equations into the mathematical model. The proposed method is developed for this specific problem and derives a set of analytical gradients to speed up convergence, and thereby having a higher computational efficiency than those typical methods.
4.3. Application in the Midcourse Guidance
In this subsection, the proposed method is used as the baseline algorithm for midcourse guidance, in which the interception scenarios are the same as that listed in Table 1. In order to demonstrate its superb performance, several typical guidance laws, which have been widely implemented for exoatmospheric interception, are also applied for comparison. Those guidance laws include proportional navigation (PN) [34], augmented proportional navigation (APN) [35], and optimal terminal guidance (OTG) [36]. It should be noted that these three methods are different from their gravity models. PN assumes that the gravity difference is zero, whereas APN and OTG assume that gravity differences are constant and linear, respectively [37].
All numerical results have been recorded, and the resulting velocity increments for various methods are listed in Table 5. Obviously, the proposed method significantly outperforms the other three guidance laws in terms of the velocity increment because it regards the minimum velocity increment as the cost function. Although the Case 1 of Group 3 is the worst case, the velocity increment of the proposed method is only 72.4% of that of PN (27.6% reduction), 73.5% of that of APN (27.5% reduction), and 73.9% of that of OTG (27.1% reduction). For some special cases, the cost reduction is even higher than 84%. Consequently, the proposed method is able to achieve very favorable fuel efficiency and retain enough energy for a successful collision. It should be noted that, because Keplerian motion is not taken into account, PN, APN, and OTG are not suitable for tail-chase interception scenarios such as the Case 3 in Group 1 and Group 2. The velocity increments for these three methods are higher than 2264 m/s and 1738 m/s for Group 1 and Group 2, respectively. However, the values of the proposed method are only 552 m/s and 421 m/s, which suggest that the proposed method has a more extensive applicability. Computational efficiency was fully discussed in Subsection 4.2. It can be concluded that the proposed method is sufficiently applicable for midcourse guidance.
5. Conclusion
In this paper, a semianalytical method is proposed for solving the exoatmospheric midcourse guidance problem with the minimum velocity increment. A nonlinear programming problem relating Lambert’s problem with Keplerian motion is firstly formulated. Then, Karush-Kuhn-Tucker conditions and its Jacobian matrix are successfully derived in an analytical manner. Therefore, Newton-Raphson method can be efficiently employed to obtain the optimal result. In order to further improve the computational efficiency, a near-optimal initialization method is developed to provide a warm start. Several experiments and the comparison with existing methods are carried out to evaluate the performance of the proposed method. Numerical results show that, even for worse case scenarios, the maximal number of iterations required is less than 6, the maximal computation time is less than 0.116 ms, and the cost reduction is more than 24.1%. Apparently, the proposed method outperforms typical methods and is of rapid convergence, high computational efficiency, and wide applicability.
Data Availability
The simulation data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.