# Optimal Midcourse Guidance Algorithm for Exoatmospheric Interception Using Analytical Gradients

Authors: Wenhao Du; Wanchun Chen; Liang Yang; Hao Zhou
Journal: International Journal of Aerospace Engineering
 (2019)

---

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


---

## Body

## 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 P-iteration method, where the slope of the TOF with respect to the semilatus rectum P 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 rectumP, is insensitive to the variation of TOF. A near-optimal initialization method, which considers P 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. Section2 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 Figure1(a). At the moment T1, 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 v0 to v1 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 Δvmin, 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).Figure 1
Illustration of the MVIP. (a) Geometry of exoatmospheric interception. (b) The MVIP is a combination of Lambert’s problem and Keplerian motion.
(a)
(b)


### 2.1. Orbital Motion of the Target

Let the orbital elements of the target orbit orb2 be(1)ELE=ΩM,iM,aM,eM,ωM,tpM,where ΩM is the right ascension of the ascending node, iM is the inclination of the orbit, aM is the semimajor axis, eM is the eccentricity, ωM is the argument of periapsis, and tpM denotes the time of periapsis passage. The subscript M 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 asfM. The target position vector, r2, and the corresponding geocentric distance, r2, are expressed as follows:
(2)r2=aM1−eM21+eMcosfM,(3)r2=r2h2=r2h2xh2yh2zT,where h2 denotes the unit vector of terminal position and is defined as follows:
(4)h2=h2xh2yh2z=cosfM+ωMcosΩM−sinfM+ωMcosiMsinΩMcosfM+ωMsinΩM+sinfM+ωMcosiMcosΩMsinfM+ωMsiniM.Furthermore, the geometric relationship betweenfM and its corresponding eccentric anomaly EM is
(5)tanEM2=1−eM1+eMtanfM2,and according to orbital motion equations, the transfer time, t, can be formulated as follows:
(6)t−tpM=aM3μEM−eMsinEM.Then, substituting equation (5) into equation (6), the orbital motion of the target can be explicitly expressed as a function of t and fM.
(7)g1t,fM=0.Equations (2) and (7) indicate that the position vector, r2, is uniquely determined by the TOF. However, equation (7) is transcendental in fM and unable to be solved analytically.


### 2.2. Lambert’s Problem

As stated above, if the terminal position vector of the target,r2, 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, r2 and v2, in terms of initial vectors, r1 and v1; thus,
(8)r2=Fr1+Gv1,v2=Ftr1+Gtv1.Clearly, the coefficients,Ft and Gt, are simply the derivatives of F and G with respect to the time, respectively, and they can be expressed in terms of cross products as follows:
(9)FGFtGt=1HEr2×v1r2×r1v2×v1v2×r1,where HE is the magnitude of HE, 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:(10)r=rcosfip+rsinfiq,(11)v=ṙcosf−rḟsinfip+ṙsinf+rḟcosfiq,where ip is a unit vector in the direction of perigee, the point of closest approach to the earth, and iq is in the orbital plane and normal to ip such that iq,ip,HE complete a right-handed coordinate system Fp. The equation of an orbit in polar coordinates is written as follows:
(12)r=P1+ecosf,where P=a1–e2 is the semilatus rectum, and a,e,f denote the semimajor axis, the eccentricity, and the true anomaly of the EKV orbit, respectively. Differentiating equation (12) and one can have
(13)ṙ=μPesinf,(14)rḟ=μP1+ecosf,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:
(15)FGFtGt=1−r2P1−cosϕr1r2sinϕμPμPtanϕ21−cosϕP−1r1−1r21−r1P1−cosϕ,where ϕ denotes the angle between r2 and r1; thus,
(16)cosϕ=r1⋅r2r1r2=r1xh2x+r1yh2y+r1zh2zr1.It should be noted that,ϕis a constant in the classical Lambert’s problem, and therefore, the Lagrangian coefficients are uniquely determined by P.In addition, the position and velocity vectors in equation (10) can be directly written in terms of the eccentric anomaly E as follows:
(17)r=acosE−eip+a1−e2sinEiq,(18)v=−aĖsinEip+a1−e2ĖcosEiq.Taking the derivative of equation (6), one can obtain that
(19)Ė=1rμa.Substituting equations (17) and (19) into equation (9), and the Lagrangian coefficients can be rewritten as follows:
(20)FGFtGt=1−ar11−cosΔEt−a3μΔE−sinΔE−μar1r2sinΔE1−ar21−cosΔE,where ΔEdenotes the difference of the eccentric anomaly, which is expressed as ΔE=E2–E1.Combining equations (15) and (20), the transfer time tcan be expressed as a function of Lagrangian coefficients
(21)t=G+a3μarccos1−r1a1−F+r1r2Ftμa,where the semimajor axis a is expressed as follows:
(22)a=μr121−F22μr11−F−Ft2r12r22.Substituting equations (22) and (15) into equation (21), Lambert’s problem can be finally expressed as a function oft and P
(23)g2t,P=0.It is apparent that equation (23) is a transcendental equation, and therefore, its solution can only be provided by iterative methods.In summary,g1 (equation (7)) and g2 (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 r2. Then, substituting the value of P into equation (15), the initial velocity vector v1 of the EKV can be calculated by Lagrangian coefficients as follows:
(24)v1=r2−Fr1G.Becausev1 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:
(25)minimizeJt=Δv2,where the square of the velocity increment is given by
(26)Δv2=v1−v02=v1x−v0x2+v1y−v0y2+v1z−v0z2.To find a local minimumt^, it is necessary to derive a set of first-order necessary conditions of this problem
(27)dJdtt=t^=∂J∂P∂P∂t+∂J∂fM∂fM∂tt=t^=0.It should be noted that two transcendental equations,g1 (equation (7)) and g2 (equation (23)), are involved in the computational process, which leads to the strong coupling among t,P,andfM 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 x=t is expanded to X=t,P,fMT. Meanwhile, g1andg2 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.


## 2.1. Orbital Motion of the Target

Let the orbital elements of the target orbit orb2 be(1)ELE=ΩM,iM,aM,eM,ωM,tpM,where ΩM is the right ascension of the ascending node, iM is the inclination of the orbit, aM is the semimajor axis, eM is the eccentricity, ωM is the argument of periapsis, and tpM denotes the time of periapsis passage. The subscript M 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 asfM. The target position vector, r2, and the corresponding geocentric distance, r2, are expressed as follows:
(2)r2=aM1−eM21+eMcosfM,(3)r2=r2h2=r2h2xh2yh2zT,where h2 denotes the unit vector of terminal position and is defined as follows:
(4)h2=h2xh2yh2z=cosfM+ωMcosΩM−sinfM+ωMcosiMsinΩMcosfM+ωMsinΩM+sinfM+ωMcosiMcosΩMsinfM+ωMsiniM.Furthermore, the geometric relationship betweenfM and its corresponding eccentric anomaly EM is
(5)tanEM2=1−eM1+eMtanfM2,and according to orbital motion equations, the transfer time, t, can be formulated as follows:
(6)t−tpM=aM3μEM−eMsinEM.Then, substituting equation (5) into equation (6), the orbital motion of the target can be explicitly expressed as a function of t and fM.
(7)g1t,fM=0.Equations (2) and (7) indicate that the position vector, r2, is uniquely determined by the TOF. However, equation (7) is transcendental in fM and unable to be solved analytically.


## 2.2. Lambert’s Problem

As stated above, if the terminal position vector of the target,r2, 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, r2 and v2, in terms of initial vectors, r1 and v1; thus,
(8)r2=Fr1+Gv1,v2=Ftr1+Gtv1.Clearly, the coefficients,Ft and Gt, are simply the derivatives of F and G with respect to the time, respectively, and they can be expressed in terms of cross products as follows:
(9)FGFtGt=1HEr2×v1r2×r1v2×v1v2×r1,where HE is the magnitude of HE, 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:(10)r=rcosfip+rsinfiq,(11)v=ṙcosf−rḟsinfip+ṙsinf+rḟcosfiq,where ip is a unit vector in the direction of perigee, the point of closest approach to the earth, and iq is in the orbital plane and normal to ip such that iq,ip,HE complete a right-handed coordinate system Fp. The equation of an orbit in polar coordinates is written as follows:
(12)r=P1+ecosf,where P=a1–e2 is the semilatus rectum, and a,e,f denote the semimajor axis, the eccentricity, and the true anomaly of the EKV orbit, respectively. Differentiating equation (12) and one can have
(13)ṙ=μPesinf,(14)rḟ=μP1+ecosf,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:
(15)FGFtGt=1−r2P1−cosϕr1r2sinϕμPμPtanϕ21−cosϕP−1r1−1r21−r1P1−cosϕ,where ϕ denotes the angle between r2 and r1; thus,
(16)cosϕ=r1⋅r2r1r2=r1xh2x+r1yh2y+r1zh2zr1.It should be noted that,ϕis a constant in the classical Lambert’s problem, and therefore, the Lagrangian coefficients are uniquely determined by P.In addition, the position and velocity vectors in equation (10) can be directly written in terms of the eccentric anomaly E as follows:
(17)r=acosE−eip+a1−e2sinEiq,(18)v=−aĖsinEip+a1−e2ĖcosEiq.Taking the derivative of equation (6), one can obtain that
(19)Ė=1rμa.Substituting equations (17) and (19) into equation (9), and the Lagrangian coefficients can be rewritten as follows:
(20)FGFtGt=1−ar11−cosΔEt−a3μΔE−sinΔE−μar1r2sinΔE1−ar21−cosΔE,where ΔEdenotes the difference of the eccentric anomaly, which is expressed as ΔE=E2–E1.Combining equations (15) and (20), the transfer time tcan be expressed as a function of Lagrangian coefficients
(21)t=G+a3μarccos1−r1a1−F+r1r2Ftμa,where the semimajor axis a is expressed as follows:
(22)a=μr121−F22μr11−F−Ft2r12r22.Substituting equations (22) and (15) into equation (21), Lambert’s problem can be finally expressed as a function oft and P
(23)g2t,P=0.It is apparent that equation (23) is a transcendental equation, and therefore, its solution can only be provided by iterative methods.In summary,g1 (equation (7)) and g2 (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 r2. Then, substituting the value of P into equation (15), the initial velocity vector v1 of the EKV can be calculated by Lagrangian coefficients as follows:
(24)v1=r2−Fr1G.Becausev1 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:
(25)minimizeJt=Δv2,where the square of the velocity increment is given by
(26)Δv2=v1−v02=v1x−v0x2+v1y−v0y2+v1z−v0z2.To find a local minimumt^, it is necessary to derive a set of first-order necessary conditions of this problem
(27)dJdtt=t^=∂J∂P∂P∂t+∂J∂fM∂fM∂tt=t^=0.It should be noted that two transcendental equations,g1 (equation (7)) and g2 (equation (23)), are involved in the computational process, which leads to the strong coupling among t,P,andfM 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 x=t is expanded to X=t,P,fMT. Meanwhile, g1andg2 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. Section3.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 functionJX subject to a set of equality constraints gX and inequality constraints mX. The problem is then
(28)minimizeJXsubjecttogX=0,mX≤0.The KKT conditions are found by defining the augmented cost function,L, so that
(29)LX,λ,η=JX−∑i=1KλigiX−∑j=1LηjmjX,where λi,ηj are the KKT multipliers associated with equality and inequality constraints, respectively. The optimal solution, X^, is determined by the gradient of the augmented cost function [33]. If J,g,andm are continuously differentiable functions, and then there exists KKT multipliers λ^ and μ^ such that
(30)∂L∂xX^,λ^,η^=0,(31)gX^=0,(32)η^mX^=0,whereη^≤0,mX^≤0.For the MVIP discussed in this paper, the corresponding NLP problem is finally expressed as follows:(33)minimizeJt,P,fM=Δv2subjecttog1t,fM=0,g2t,P,fM=0.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 ing1. Thus, fM is involved in the third equation of equation (33), which also decreases the computational difficulty. Therefore, the independent variables in the equality constraint g2 are set to be X=t,P,fMT, which is different from that in equation (23). Furthermore, the augmented cost function is given by
(34)L=J−λ1g1−λ2g2.Taking the partial derivative of the independent variables, one can obtain(35)∂L∂X=0,X=t,P,fM.Rearranging equation (35), then the KKT multipliers λ1 and λ2 can be formulated as follows:
(36)λ1=−aM3μλ2,(37)λ2=∂J/∂P∂g2/∂P.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:
(38)ΦX=g1g2g3T,where
(39)g1t,fM=EM−eMsinEM−μaM3t−tpM,(40)g2t,P,fM=G+a3μarccos1−r1a1−F+r1r2Ftμa−t,(41)g3P,fM=∂J∂fM+aM3μ∂J/∂P∂g2/∂P∂g1∂fM−∂J/∂P∂g2/∂P∂g2∂fM,where g1 and g2 describe the mathematical models of Keplerian motion and Lambert’s problem in the MVIP. Meanwhile, g3 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:(42)XK+1=XK−ΦXXK−1ΦXK,where the superscript K means the K-th iteration, and the Jacobian matrix ΦXX is given by
(43)ΦXX=∂ΦX∂X=∂g1∂t∂g1∂P∂g1∂fM∂g2∂t∂g2∂P∂g2∂fM∂g3∂t∂g3∂P∂g3∂fM.In spite of mathematical complication, analytical expressions of all the parameters in theΦXX 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 ofg1 are firstly formulated due to its concise structure.


#### 3.2.1. Partial Derivatives ofg1

SinceP is not explicit in g1, the partial derivative of g1 with respect to P is zero.
(44)∂g1∂P=0.The partial derivative ofg1 with respect to tcan be obviously expressed as follows:
(45)∂g1∂t=−μaM3.Substituting equation (5) into equation (39), the partial derivative of g1 with respect to fM is given according to trigonometry formulas
(46)∂g1∂fM=1−eM23/21+eMcosfM2.Although the partial derivatives ofg1 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 ofg2

Observing the expression ofg2 in equation (39), the partial derivative with respect to t can be obviously expressed as follows:
(47)∂g2∂t=−1.For simplifying the derivations, three notations are introduced as follows:(48)k=r1r21−cosϕ,l=r1+r2,m=r1r21+cosϕ.Substituting equation (15) into equation (22), the semimajor axis of the intercept orbit is rewritten as follows:
(49)a=mkP2m−l2P2+2klP−k2.Then, the partial derivative ofawith respect to P is expressed as follows:
(50)∂a∂P=mkl2−2mP2−k2χ02,where χ0 is another notation expressed as follows:
(51)χ0=2m−l2P2+2klP−k2.Substituting equations (15) and (49) into equation (39) and partially differentiating g2 with respect to P yields
(52)∂g2∂P=−G2P+32aμ∂a∂Pχ3+a3μχ4,where
(53)χ1=1−kaP,χ2=2aPk−k2aP,χ3=arccosχ1+tanϕ/2k/P−laP,χ4=2kχ2+aPkP−ltanϕ2k−PlmP3−kP2aPtanϕ2.Partially differentiating equation (2) with respect to fM, one can obtain that
(54)∂r2∂fM=r2eMsinfM1+eMcosfM.Substituting equation (16) into the expression of G in equation (15), it can be rewritten as follows:
(55)G=r2μPr12−r1xh2x+r1yh2y+r1zh2z2.Partially differentiating equation (55) with respect tofM, and introducing notations as follows:
(56)r1h2=r1xh2x+r1yh2y+r1zh2z,r1ph2=r1x∂h2x∂fM+r1y∂h2y∂fM+r1z∂h2z∂fM,sr12h2=r12−r1xh2x+r1yh2y+r1zh2z2,where
(57)∂h2∂fM=∂h2x∂fM∂h2y∂fM∂h2z∂fM=−sinfM+wcosΩ−cosfM+wcosisinΩ−sinfM+wsinΩ+cosfM+wcosicosΩcosfM+wsini.Then, the partial derivative of the Lagrangian coefficientG with respect to fM is given by
(58)∂G∂fM=1μP∂r2∂fMsr12h2−r2r1h2r1ph2sr12h2.In order to simplify the expression, notationsχ5,χ51,χ6,χ61~χ65,χ7, and χ8 are also introduced, which are expressed as follows:
(59)χ5=2r2∂r2∂fMsr12h22−r2r1h2r1ph2,χ51=2P2χ61+2Pχ62−2r2χ6χ63,χ6=r1−r1h2,χ61=r1h2−r2∂r2∂fM+r2r1ph2,χ62=r1+2r2∂r2∂fMχ6−r1+r2r2r1ph2,χ63=χ6∂r2∂fM−r2r1ph2,χ64=mkχ02,χ65=χ63−r2χ6a∂a∂fM,χ7=χ6sr12h2,χ8=r2χ6P−r2−r1.Partially differentiating equation (49) with respect to fM, and substituting equation (54) into the result, one finds
(60)∂a∂fM=Pχ5χ0−χ64⋅χ51.Then, substituting equations (15) and (55) into equation (39) and partially differentiating g2 with respect to fM yields
(61)∂g2∂fM=∂G∂fM+32aμ∂a∂fMarccosχ1+aμχ65Pχ2+1μPχ7χ8∂a∂fM+aχ8∂χ7∂fM+aχ7∂χ8∂fM,where
(62)∂χ7∂fM=−r1r1ph2r1+r1h2sr12h2,∂χ8∂fM=χ63P−∂r2∂fM.So far, the partial derivatives ofg2 have been successfully expressed by equations (47), (52), and (61) in an analytical manner. It is worth noting that g3 consists of the partial derivatives of J,g1,andg2. Therefore, the partial derivatives of g3 involve a number of second-order partial derivatives, which are much more complicated.


#### 3.2.3. Partial Derivatives ofg3

Observing the expression forg3 in equation (39), it is easy to see that t does not appear explicitly, so the partial derivative of g3 with respect to t can be expressed as follows:
(63)∂g3∂t=0.Then, let us focus on the partial derivative with respect toP. The derivation is presented as follows:At the beginning, substituting equations (2) and (24) into equation (33) and introducing a set of notations
(64)JPx1=r2xG−FGr1x−v0x,JPy1=r2yG−FGr1y−v0y,JPz1=r2zG−FGr1z−v0z,then the cost function can be rewritten as follows:
(65)J=JPx12+JPy12+JPz12.Partially differentiating equation (65) with respect to P and introducing another set of notations as
(66)JPx2=r2x−r1x2−F2PG,JPy2=r2y−r1y2−F2PG,JPz2=r2z−r1z2−F2PG,one has
(67)∂J∂P=2JPx1JPx2+JPy1JPy2+JPz1JPz2.Then, differentiating equations (50), (67), and (52) with respect toP, one can obtain second-order partial derivatives
(68)∂2a∂P2=2mk2m−l22P3+32m−l2Pk2+2k3lχ03,∂2J∂P2=2JPx22+JPy22+JPz22+1−FP2Gr1xJPx1+r1yJPy1+r1zJPz1−12P∂J∂P,∂2g2∂P2=34GP2+32121aμ∂a∂P2+aμ∂2a∂P2χ3+3aμ∂a∂Pχ4+a3μ∂χ4∂P,where
(69)∂χ4∂P=2kχ2+aP2Pl−3kmP4−4ak−Pl2kχ1m2P5χ2+tanϕ/2mP4aP2kak−aPl−mP+k−Pl22P∂a∂P+a.Partially differentiatingF in equation (15) and using equations (16) and (54), one has
(70)∂F∂fM=cosϕ−1P∂r2∂fM+r2⋅r1ph2Pr1.Then, differentiating equations (50), (65), and (52) with respect to fM and introducing two sets of notations expressed as
(71)JfMx1=∂r2∂fMh2x+r2∂h2x∂fMg−r2h2x∂G∂fM,JfMx2=h2xG∂r2∂fM+r2G∂h2x∂fM−r2h2x−r1xFG2∂G∂fM−r1xG∂F∂fM,JfMy1=∂r2∂fMh2y+r2∂h2y∂fMg−r2h2y∂G∂fM,JfMy2=h2yG∂r2∂fM+r2G∂h2y∂fM−r2h2y−r1yFG2∂G∂fM−r1yG∂F∂fM,JfMz1=∂r2∂fMh2z+r2∂h2z∂fMg−r2h2z∂G∂fM,JfMz2=h2zG∂r2∂fM+r2G∂h2z∂fM−r2h2z−r1zFG2∂G∂fM−r1zG∂F∂fM,one can obtain the second-order partial and mixed derivatives
(72)∂2a∂P∂fM=∂a∂fM1P−2Pχ032m−l2P+klχ0χ5−2mkχ51+2Pχ61+χ62mkχ0,∂2J∂P∂fM=2JPx12PG2JfMx1+Gr1x∂F∂fM+2−Fr1x∂G∂fM+JPx2JfMx2+JPy12PG2JfMy1+Gr1y∂F∂fM+2−Fr1y∂G∂fM+JPy2JfMy2+JPz12PG2JfMz1+Gr1z∂F∂fM+2−Fr1z∂G∂fM+JPz2JfMz2,∂2g2∂P∂fM=−12P∂G∂fM+32arccosχ1−r2χ6aPχ2aμ∂2a∂P∂fM+k−Pl3k−2r2χ6mP3aμ∂a∂fM−aμχ63P2χ21+χ1+12aμ∂a∂P32∂a∂fMarccosχ1+χ65Pχ2−r2χ6χ1aμP2χ21+χ1.Substituting equations (46), (52), (61), and (67) into equation (39), and partially differentiating with respect to P, the partial derivative of g3 with respect to P is finally derived as follows:
(73)∂g3∂P=∂2J∂P∂fM+∂2J∂P2∂g2∂P−∂J∂P∂2g2∂P2⋅aM3/μ∂g1/∂fM−∂g2/∂fM∂g2/∂P2−∂J/∂P∂g2/∂g2⋅∂2g2∂P∂fM.Let us focus on the last term of the partial derivative ofg3. Partially differentiating equations (54), (58) ,and (70) with respect to fM, one can obtain that
(74)∂2r2∂fM2=2e2−e2cos2fM+ecosfM1+ecosfM2r2,∂2G∂fM2=sr12h2μP∂2r2∂fM2−2r1ph2r1h2sr12h22∂r2∂fM−r1ph22−r1h22r12+r1h24sr12h23r2,∂2F∂fM2=1r1P2∂r2∂fMr1ph2−r2r1h2−∂2r2∂fM2r1−r1h2.Partially differentiating equation (65) with respect to fM twice and introducing another set of notations are expressed as follows:
(75)JfMx3=∂2r2∂fM2h2x+2∂r2∂fM∂h2x∂fM−r2h2x1G−r2h2xG2∂2G∂fM2−2JfMx1G∂G∂fM−r1xJfMFG,JfMy3=∂2r2∂fM2⋅h2y+2∂r2∂fM∂h2y∂fM−r2h2y1G−r2h2yG2∂2G∂fM2−2JfMy1G∂G∂fM−r1yJfMFG,JfMz3=∂2r2∂fM2⋅h2z+2∂r2∂fM∂h2z∂fM−r2h2z1G−r2h2zG2∂2G∂fM2−2JfMz1G∂G∂fM−r1zJfMFG,where
(76)JfMFG=∂2F∂fM21G−FG2∂2G∂fM2−2G2∂F∂fM∂G∂fM+2FG3∂G∂fM2.Then, the second-order partial derivative of the cost functionJ with respect to fM is given by
(77)∂2J∂fM2=2JfMx22+JPx1JfMx3+JfMy22+JPy1JfMy3+JfMz22+JPz1JfMz3.Partially differentiating equation (60) with respect to fM, then one has
(78)∂2a∂fM2=P∂χ5∂fM1χ0−χ64∂χ51∂fM−2χ5χ51χ02+2χ64χ512χ0,where
(79)χ9=∂2r2∂fM2χ6−2∂r2∂fMr1ph2+r2r1h2,∂χ5∂fM=∂r2∂fM2+r2∂2r2∂fM22sr12h22−8r2r1h2r1ph2∂r2∂fM+2r22r1h22−r1ph22,∂χ51∂fM=r1−r2∂2r2∂fM2−χ9−∂r2∂fM22P2−2χ632−2r2χ6χ9+r1+2r2χ9−r2r1h2+2∂r2∂fM2χ62P.Partially differentiating equation (61), one can obtain the second-order partial derivative of g2 with respect to fM
(80)∂2g2∂fM2=∂2G∂fM2+32aμ∂2a∂fM2+12a∂a∂fM2arccosχ1−χ65Pχ2aμ∂a∂fM+aμχ91+1μP∂2a∂fM2χ7χ8+a∂2χ7∂fM2χ8+aχ7∂2χ8∂fM2+2∂a∂fM∂χ7∂fMχ8+2∂a∂fMχ7∂χ8∂fM+2a∂χ7∂fM∂χ8∂fM,where
(81)∂2χ7∂fM2=−r1h2r1ph2+r1ph2sr12h2r1h2sr12h2−χ7∂χ7∂fM,∂2χ8∂fM2=χ9P−∂2r2∂fM2.Furthermore, partially differentiating equation (46), the second-order partial derivative of g1 is expressed as follows:
(82)∂2g1∂fM2=eM1−eM2sinfM1+eMcosfM+eMsinEM.Then, the partial derivative ofg3 with respect to fM can be expressed as follows:
(83)∂g3∂fM=∂2J∂fM2+∂2J∂P∂fM∂g2∂P−∂J∂P∂2g2∂P∂fMaM3/μ∂g1/∂fM−∂g2/∂fM∂g2/∂P2+∂J/∂P∂g2/∂PaM3μ∂2g1∂fM2−∂2g2∂fM2.Now, the partial derivatives ofg3 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 variablesX=t,P,fMT for the MVIP. This method is based on a discovery that semilatus rectum, P, is insensitive to the variation of TOF in the computation process. In fact, the semilatus rectum P 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 P 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,Pis 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, t, and the true anomaly, fM. It should be noted that, according to equation (7), t can be analytically expressed as a function of fM. Therefore, the dimension of the search space can be further reduced to formulate a single-variable rooting problem. Let ΨfM denote the equality constraint on fM, which is expressed as follows:
(84)ΨfM=g2−aM3μg1.To determine the initial guess offM, a single-variable Newton-Raphson iteration equation is formulated as follows:
(85)fML+1=fML−ΨfMLΨ̇fML,where
(86)Ψ̇fM=dg2dfM−aM3μdg1dfM.It should be pointed out thatΨ̇fM can be calculated analytically from the derivation in Section 3.2.The iteration is finished if the difference betweentL and tL+1 is small enough, and the convergence criterion of the initialization method is given by
(87)tL+1−tL<εt.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εX, 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:
(88)XK+1−XK<εX.Otherwise, the iterative variablesX=t,P,fMT 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 rectumP reaches as high as 106, whereas the true anomaly fM 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 u and the corresponding convergence criterion is formulated as follows:
(89)uXK+1−uXK<εx,where u is the weight matrix.The procedure of implementing the proposed algorithm is included in the flowchart in Figure2. 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.Figure 2
Flowchart of the proposed method.


## 3.1. Derivation of the KKT Conditions

Consider a NLP problem with both equality and inequality constraints. The objective is to minimize the cost functionJX subject to a set of equality constraints gX and inequality constraints mX. The problem is then
(28)minimizeJXsubjecttogX=0,mX≤0.The KKT conditions are found by defining the augmented cost function,L, so that
(29)LX,λ,η=JX−∑i=1KλigiX−∑j=1LηjmjX,where λi,ηj are the KKT multipliers associated with equality and inequality constraints, respectively. The optimal solution, X^, is determined by the gradient of the augmented cost function [33]. If J,g,andm are continuously differentiable functions, and then there exists KKT multipliers λ^ and μ^ such that
(30)∂L∂xX^,λ^,η^=0,(31)gX^=0,(32)η^mX^=0,whereη^≤0,mX^≤0.For the MVIP discussed in this paper, the corresponding NLP problem is finally expressed as follows:(33)minimizeJt,P,fM=Δv2subjecttog1t,fM=0,g2t,P,fM=0.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 ing1. Thus, fM is involved in the third equation of equation (33), which also decreases the computational difficulty. Therefore, the independent variables in the equality constraint g2 are set to be X=t,P,fMT, which is different from that in equation (23). Furthermore, the augmented cost function is given by
(34)L=J−λ1g1−λ2g2.Taking the partial derivative of the independent variables, one can obtain(35)∂L∂X=0,X=t,P,fM.Rearranging equation (35), then the KKT multipliers λ1 and λ2 can be formulated as follows:
(36)λ1=−aM3μλ2,(37)λ2=∂J/∂P∂g2/∂P.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:
(38)ΦX=g1g2g3T,where
(39)g1t,fM=EM−eMsinEM−μaM3t−tpM,(40)g2t,P,fM=G+a3μarccos1−r1a1−F+r1r2Ftμa−t,(41)g3P,fM=∂J∂fM+aM3μ∂J/∂P∂g2/∂P∂g1∂fM−∂J/∂P∂g2/∂P∂g2∂fM,where g1 and g2 describe the mathematical models of Keplerian motion and Lambert’s problem in the MVIP. Meanwhile, g3 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:(42)XK+1=XK−ΦXXK−1ΦXK,where the superscript K means the K-th iteration, and the Jacobian matrix ΦXX is given by
(43)ΦXX=∂ΦX∂X=∂g1∂t∂g1∂P∂g1∂fM∂g2∂t∂g2∂P∂g2∂fM∂g3∂t∂g3∂P∂g3∂fM.In spite of mathematical complication, analytical expressions of all the parameters in theΦXX 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 ofg1 are firstly formulated due to its concise structure.


### 3.2.1. Partial Derivatives ofg1

SinceP is not explicit in g1, the partial derivative of g1 with respect to P is zero.
(44)∂g1∂P=0.The partial derivative ofg1 with respect to tcan be obviously expressed as follows:
(45)∂g1∂t=−μaM3.Substituting equation (5) into equation (39), the partial derivative of g1 with respect to fM is given according to trigonometry formulas
(46)∂g1∂fM=1−eM23/21+eMcosfM2.Although the partial derivatives ofg1 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 ofg2

Observing the expression ofg2 in equation (39), the partial derivative with respect to t can be obviously expressed as follows:
(47)∂g2∂t=−1.For simplifying the derivations, three notations are introduced as follows:(48)k=r1r21−cosϕ,l=r1+r2,m=r1r21+cosϕ.Substituting equation (15) into equation (22), the semimajor axis of the intercept orbit is rewritten as follows:
(49)a=mkP2m−l2P2+2klP−k2.Then, the partial derivative ofawith respect to P is expressed as follows:
(50)∂a∂P=mkl2−2mP2−k2χ02,where χ0 is another notation expressed as follows:
(51)χ0=2m−l2P2+2klP−k2.Substituting equations (15) and (49) into equation (39) and partially differentiating g2 with respect to P yields
(52)∂g2∂P=−G2P+32aμ∂a∂Pχ3+a3μχ4,where
(53)χ1=1−kaP,χ2=2aPk−k2aP,χ3=arccosχ1+tanϕ/2k/P−laP,χ4=2kχ2+aPkP−ltanϕ2k−PlmP3−kP2aPtanϕ2.Partially differentiating equation (2) with respect to fM, one can obtain that
(54)∂r2∂fM=r2eMsinfM1+eMcosfM.Substituting equation (16) into the expression of G in equation (15), it can be rewritten as follows:
(55)G=r2μPr12−r1xh2x+r1yh2y+r1zh2z2.Partially differentiating equation (55) with respect tofM, and introducing notations as follows:
(56)r1h2=r1xh2x+r1yh2y+r1zh2z,r1ph2=r1x∂h2x∂fM+r1y∂h2y∂fM+r1z∂h2z∂fM,sr12h2=r12−r1xh2x+r1yh2y+r1zh2z2,where
(57)∂h2∂fM=∂h2x∂fM∂h2y∂fM∂h2z∂fM=−sinfM+wcosΩ−cosfM+wcosisinΩ−sinfM+wsinΩ+cosfM+wcosicosΩcosfM+wsini.Then, the partial derivative of the Lagrangian coefficientG with respect to fM is given by
(58)∂G∂fM=1μP∂r2∂fMsr12h2−r2r1h2r1ph2sr12h2.In order to simplify the expression, notationsχ5,χ51,χ6,χ61~χ65,χ7, and χ8 are also introduced, which are expressed as follows:
(59)χ5=2r2∂r2∂fMsr12h22−r2r1h2r1ph2,χ51=2P2χ61+2Pχ62−2r2χ6χ63,χ6=r1−r1h2,χ61=r1h2−r2∂r2∂fM+r2r1ph2,χ62=r1+2r2∂r2∂fMχ6−r1+r2r2r1ph2,χ63=χ6∂r2∂fM−r2r1ph2,χ64=mkχ02,χ65=χ63−r2χ6a∂a∂fM,χ7=χ6sr12h2,χ8=r2χ6P−r2−r1.Partially differentiating equation (49) with respect to fM, and substituting equation (54) into the result, one finds
(60)∂a∂fM=Pχ5χ0−χ64⋅χ51.Then, substituting equations (15) and (55) into equation (39) and partially differentiating g2 with respect to fM yields
(61)∂g2∂fM=∂G∂fM+32aμ∂a∂fMarccosχ1+aμχ65Pχ2+1μPχ7χ8∂a∂fM+aχ8∂χ7∂fM+aχ7∂χ8∂fM,where
(62)∂χ7∂fM=−r1r1ph2r1+r1h2sr12h2,∂χ8∂fM=χ63P−∂r2∂fM.So far, the partial derivatives ofg2 have been successfully expressed by equations (47), (52), and (61) in an analytical manner. It is worth noting that g3 consists of the partial derivatives of J,g1,andg2. Therefore, the partial derivatives of g3 involve a number of second-order partial derivatives, which are much more complicated.


### 3.2.3. Partial Derivatives ofg3

Observing the expression forg3 in equation (39), it is easy to see that t does not appear explicitly, so the partial derivative of g3 with respect to t can be expressed as follows:
(63)∂g3∂t=0.Then, let us focus on the partial derivative with respect toP. The derivation is presented as follows:At the beginning, substituting equations (2) and (24) into equation (33) and introducing a set of notations
(64)JPx1=r2xG−FGr1x−v0x,JPy1=r2yG−FGr1y−v0y,JPz1=r2zG−FGr1z−v0z,then the cost function can be rewritten as follows:
(65)J=JPx12+JPy12+JPz12.Partially differentiating equation (65) with respect to P and introducing another set of notations as
(66)JPx2=r2x−r1x2−F2PG,JPy2=r2y−r1y2−F2PG,JPz2=r2z−r1z2−F2PG,one has
(67)∂J∂P=2JPx1JPx2+JPy1JPy2+JPz1JPz2.Then, differentiating equations (50), (67), and (52) with respect toP, one can obtain second-order partial derivatives
(68)∂2a∂P2=2mk2m−l22P3+32m−l2Pk2+2k3lχ03,∂2J∂P2=2JPx22+JPy22+JPz22+1−FP2Gr1xJPx1+r1yJPy1+r1zJPz1−12P∂J∂P,∂2g2∂P2=34GP2+32121aμ∂a∂P2+aμ∂2a∂P2χ3+3aμ∂a∂Pχ4+a3μ∂χ4∂P,where
(69)∂χ4∂P=2kχ2+aP2Pl−3kmP4−4ak−Pl2kχ1m2P5χ2+tanϕ/2mP4aP2kak−aPl−mP+k−Pl22P∂a∂P+a.Partially differentiatingF in equation (15) and using equations (16) and (54), one has
(70)∂F∂fM=cosϕ−1P∂r2∂fM+r2⋅r1ph2Pr1.Then, differentiating equations (50), (65), and (52) with respect to fM and introducing two sets of notations expressed as
(71)JfMx1=∂r2∂fMh2x+r2∂h2x∂fMg−r2h2x∂G∂fM,JfMx2=h2xG∂r2∂fM+r2G∂h2x∂fM−r2h2x−r1xFG2∂G∂fM−r1xG∂F∂fM,JfMy1=∂r2∂fMh2y+r2∂h2y∂fMg−r2h2y∂G∂fM,JfMy2=h2yG∂r2∂fM+r2G∂h2y∂fM−r2h2y−r1yFG2∂G∂fM−r1yG∂F∂fM,JfMz1=∂r2∂fMh2z+r2∂h2z∂fMg−r2h2z∂G∂fM,JfMz2=h2zG∂r2∂fM+r2G∂h2z∂fM−r2h2z−r1zFG2∂G∂fM−r1zG∂F∂fM,one can obtain the second-order partial and mixed derivatives
(72)∂2a∂P∂fM=∂a∂fM1P−2Pχ032m−l2P+klχ0χ5−2mkχ51+2Pχ61+χ62mkχ0,∂2J∂P∂fM=2JPx12PG2JfMx1+Gr1x∂F∂fM+2−Fr1x∂G∂fM+JPx2JfMx2+JPy12PG2JfMy1+Gr1y∂F∂fM+2−Fr1y∂G∂fM+JPy2JfMy2+JPz12PG2JfMz1+Gr1z∂F∂fM+2−Fr1z∂G∂fM+JPz2JfMz2,∂2g2∂P∂fM=−12P∂G∂fM+32arccosχ1−r2χ6aPχ2aμ∂2a∂P∂fM+k−Pl3k−2r2χ6mP3aμ∂a∂fM−aμχ63P2χ21+χ1+12aμ∂a∂P32∂a∂fMarccosχ1+χ65Pχ2−r2χ6χ1aμP2χ21+χ1.Substituting equations (46), (52), (61), and (67) into equation (39), and partially differentiating with respect to P, the partial derivative of g3 with respect to P is finally derived as follows:
(73)∂g3∂P=∂2J∂P∂fM+∂2J∂P2∂g2∂P−∂J∂P∂2g2∂P2⋅aM3/μ∂g1/∂fM−∂g2/∂fM∂g2/∂P2−∂J/∂P∂g2/∂g2⋅∂2g2∂P∂fM.Let us focus on the last term of the partial derivative ofg3. Partially differentiating equations (54), (58) ,and (70) with respect to fM, one can obtain that
(74)∂2r2∂fM2=2e2−e2cos2fM+ecosfM1+ecosfM2r2,∂2G∂fM2=sr12h2μP∂2r2∂fM2−2r1ph2r1h2sr12h22∂r2∂fM−r1ph22−r1h22r12+r1h24sr12h23r2,∂2F∂fM2=1r1P2∂r2∂fMr1ph2−r2r1h2−∂2r2∂fM2r1−r1h2.Partially differentiating equation (65) with respect to fM twice and introducing another set of notations are expressed as follows:
(75)JfMx3=∂2r2∂fM2h2x+2∂r2∂fM∂h2x∂fM−r2h2x1G−r2h2xG2∂2G∂fM2−2JfMx1G∂G∂fM−r1xJfMFG,JfMy3=∂2r2∂fM2⋅h2y+2∂r2∂fM∂h2y∂fM−r2h2y1G−r2h2yG2∂2G∂fM2−2JfMy1G∂G∂fM−r1yJfMFG,JfMz3=∂2r2∂fM2⋅h2z+2∂r2∂fM∂h2z∂fM−r2h2z1G−r2h2zG2∂2G∂fM2−2JfMz1G∂G∂fM−r1zJfMFG,where
(76)JfMFG=∂2F∂fM21G−FG2∂2G∂fM2−2G2∂F∂fM∂G∂fM+2FG3∂G∂fM2.Then, the second-order partial derivative of the cost functionJ with respect to fM is given by
(77)∂2J∂fM2=2JfMx22+JPx1JfMx3+JfMy22+JPy1JfMy3+JfMz22+JPz1JfMz3.Partially differentiating equation (60) with respect to fM, then one has
(78)∂2a∂fM2=P∂χ5∂fM1χ0−χ64∂χ51∂fM−2χ5χ51χ02+2χ64χ512χ0,where
(79)χ9=∂2r2∂fM2χ6−2∂r2∂fMr1ph2+r2r1h2,∂χ5∂fM=∂r2∂fM2+r2∂2r2∂fM22sr12h22−8r2r1h2r1ph2∂r2∂fM+2r22r1h22−r1ph22,∂χ51∂fM=r1−r2∂2r2∂fM2−χ9−∂r2∂fM22P2−2χ632−2r2χ6χ9+r1+2r2χ9−r2r1h2+2∂r2∂fM2χ62P.Partially differentiating equation (61), one can obtain the second-order partial derivative of g2 with respect to fM
(80)∂2g2∂fM2=∂2G∂fM2+32aμ∂2a∂fM2+12a∂a∂fM2arccosχ1−χ65Pχ2aμ∂a∂fM+aμχ91+1μP∂2a∂fM2χ7χ8+a∂2χ7∂fM2χ8+aχ7∂2χ8∂fM2+2∂a∂fM∂χ7∂fMχ8+2∂a∂fMχ7∂χ8∂fM+2a∂χ7∂fM∂χ8∂fM,where
(81)∂2χ7∂fM2=−r1h2r1ph2+r1ph2sr12h2r1h2sr12h2−χ7∂χ7∂fM,∂2χ8∂fM2=χ9P−∂2r2∂fM2.Furthermore, partially differentiating equation (46), the second-order partial derivative of g1 is expressed as follows:
(82)∂2g1∂fM2=eM1−eM2sinfM1+eMcosfM+eMsinEM.Then, the partial derivative ofg3 with respect to fM can be expressed as follows:
(83)∂g3∂fM=∂2J∂fM2+∂2J∂P∂fM∂g2∂P−∂J∂P∂2g2∂P∂fMaM3/μ∂g1/∂fM−∂g2/∂fM∂g2/∂P2+∂J/∂P∂g2/∂PaM3μ∂2g1∂fM2−∂2g2∂fM2.Now, the partial derivatives ofg3 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.2.1. Partial Derivatives ofg1

SinceP is not explicit in g1, the partial derivative of g1 with respect to P is zero.
(44)∂g1∂P=0.The partial derivative ofg1 with respect to tcan be obviously expressed as follows:
(45)∂g1∂t=−μaM3.Substituting equation (5) into equation (39), the partial derivative of g1 with respect to fM is given according to trigonometry formulas
(46)∂g1∂fM=1−eM23/21+eMcosfM2.Although the partial derivatives ofg1 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 ofg2

Observing the expression ofg2 in equation (39), the partial derivative with respect to t can be obviously expressed as follows:
(47)∂g2∂t=−1.For simplifying the derivations, three notations are introduced as follows:(48)k=r1r21−cosϕ,l=r1+r2,m=r1r21+cosϕ.Substituting equation (15) into equation (22), the semimajor axis of the intercept orbit is rewritten as follows:
(49)a=mkP2m−l2P2+2klP−k2.Then, the partial derivative ofawith respect to P is expressed as follows:
(50)∂a∂P=mkl2−2mP2−k2χ02,where χ0 is another notation expressed as follows:
(51)χ0=2m−l2P2+2klP−k2.Substituting equations (15) and (49) into equation (39) and partially differentiating g2 with respect to P yields
(52)∂g2∂P=−G2P+32aμ∂a∂Pχ3+a3μχ4,where
(53)χ1=1−kaP,χ2=2aPk−k2aP,χ3=arccosχ1+tanϕ/2k/P−laP,χ4=2kχ2+aPkP−ltanϕ2k−PlmP3−kP2aPtanϕ2.Partially differentiating equation (2) with respect to fM, one can obtain that
(54)∂r2∂fM=r2eMsinfM1+eMcosfM.Substituting equation (16) into the expression of G in equation (15), it can be rewritten as follows:
(55)G=r2μPr12−r1xh2x+r1yh2y+r1zh2z2.Partially differentiating equation (55) with respect tofM, and introducing notations as follows:
(56)r1h2=r1xh2x+r1yh2y+r1zh2z,r1ph2=r1x∂h2x∂fM+r1y∂h2y∂fM+r1z∂h2z∂fM,sr12h2=r12−r1xh2x+r1yh2y+r1zh2z2,where
(57)∂h2∂fM=∂h2x∂fM∂h2y∂fM∂h2z∂fM=−sinfM+wcosΩ−cosfM+wcosisinΩ−sinfM+wsinΩ+cosfM+wcosicosΩcosfM+wsini.Then, the partial derivative of the Lagrangian coefficientG with respect to fM is given by
(58)∂G∂fM=1μP∂r2∂fMsr12h2−r2r1h2r1ph2sr12h2.In order to simplify the expression, notationsχ5,χ51,χ6,χ61~χ65,χ7, and χ8 are also introduced, which are expressed as follows:
(59)χ5=2r2∂r2∂fMsr12h22−r2r1h2r1ph2,χ51=2P2χ61+2Pχ62−2r2χ6χ63,χ6=r1−r1h2,χ61=r1h2−r2∂r2∂fM+r2r1ph2,χ62=r1+2r2∂r2∂fMχ6−r1+r2r2r1ph2,χ63=χ6∂r2∂fM−r2r1ph2,χ64=mkχ02,χ65=χ63−r2χ6a∂a∂fM,χ7=χ6sr12h2,χ8=r2χ6P−r2−r1.Partially differentiating equation (49) with respect to fM, and substituting equation (54) into the result, one finds
(60)∂a∂fM=Pχ5χ0−χ64⋅χ51.Then, substituting equations (15) and (55) into equation (39) and partially differentiating g2 with respect to fM yields
(61)∂g2∂fM=∂G∂fM+32aμ∂a∂fMarccosχ1+aμχ65Pχ2+1μPχ7χ8∂a∂fM+aχ8∂χ7∂fM+aχ7∂χ8∂fM,where
(62)∂χ7∂fM=−r1r1ph2r1+r1h2sr12h2,∂χ8∂fM=χ63P−∂r2∂fM.So far, the partial derivatives ofg2 have been successfully expressed by equations (47), (52), and (61) in an analytical manner. It is worth noting that g3 consists of the partial derivatives of J,g1,andg2. Therefore, the partial derivatives of g3 involve a number of second-order partial derivatives, which are much more complicated.


## 3.2.3. Partial Derivatives ofg3

Observing the expression forg3 in equation (39), it is easy to see that t does not appear explicitly, so the partial derivative of g3 with respect to t can be expressed as follows:
(63)∂g3∂t=0.Then, let us focus on the partial derivative with respect toP. The derivation is presented as follows:At the beginning, substituting equations (2) and (24) into equation (33) and introducing a set of notations
(64)JPx1=r2xG−FGr1x−v0x,JPy1=r2yG−FGr1y−v0y,JPz1=r2zG−FGr1z−v0z,then the cost function can be rewritten as follows:
(65)J=JPx12+JPy12+JPz12.Partially differentiating equation (65) with respect to P and introducing another set of notations as
(66)JPx2=r2x−r1x2−F2PG,JPy2=r2y−r1y2−F2PG,JPz2=r2z−r1z2−F2PG,one has
(67)∂J∂P=2JPx1JPx2+JPy1JPy2+JPz1JPz2.Then, differentiating equations (50), (67), and (52) with respect toP, one can obtain second-order partial derivatives
(68)∂2a∂P2=2mk2m−l22P3+32m−l2Pk2+2k3lχ03,∂2J∂P2=2JPx22+JPy22+JPz22+1−FP2Gr1xJPx1+r1yJPy1+r1zJPz1−12P∂J∂P,∂2g2∂P2=34GP2+32121aμ∂a∂P2+aμ∂2a∂P2χ3+3aμ∂a∂Pχ4+a3μ∂χ4∂P,where
(69)∂χ4∂P=2kχ2+aP2Pl−3kmP4−4ak−Pl2kχ1m2P5χ2+tanϕ/2mP4aP2kak−aPl−mP+k−Pl22P∂a∂P+a.Partially differentiatingF in equation (15) and using equations (16) and (54), one has
(70)∂F∂fM=cosϕ−1P∂r2∂fM+r2⋅r1ph2Pr1.Then, differentiating equations (50), (65), and (52) with respect to fM and introducing two sets of notations expressed as
(71)JfMx1=∂r2∂fMh2x+r2∂h2x∂fMg−r2h2x∂G∂fM,JfMx2=h2xG∂r2∂fM+r2G∂h2x∂fM−r2h2x−r1xFG2∂G∂fM−r1xG∂F∂fM,JfMy1=∂r2∂fMh2y+r2∂h2y∂fMg−r2h2y∂G∂fM,JfMy2=h2yG∂r2∂fM+r2G∂h2y∂fM−r2h2y−r1yFG2∂G∂fM−r1yG∂F∂fM,JfMz1=∂r2∂fMh2z+r2∂h2z∂fMg−r2h2z∂G∂fM,JfMz2=h2zG∂r2∂fM+r2G∂h2z∂fM−r2h2z−r1zFG2∂G∂fM−r1zG∂F∂fM,one can obtain the second-order partial and mixed derivatives
(72)∂2a∂P∂fM=∂a∂fM1P−2Pχ032m−l2P+klχ0χ5−2mkχ51+2Pχ61+χ62mkχ0,∂2J∂P∂fM=2JPx12PG2JfMx1+Gr1x∂F∂fM+2−Fr1x∂G∂fM+JPx2JfMx2+JPy12PG2JfMy1+Gr1y∂F∂fM+2−Fr1y∂G∂fM+JPy2JfMy2+JPz12PG2JfMz1+Gr1z∂F∂fM+2−Fr1z∂G∂fM+JPz2JfMz2,∂2g2∂P∂fM=−12P∂G∂fM+32arccosχ1−r2χ6aPχ2aμ∂2a∂P∂fM+k−Pl3k−2r2χ6mP3aμ∂a∂fM−aμχ63P2χ21+χ1+12aμ∂a∂P32∂a∂fMarccosχ1+χ65Pχ2−r2χ6χ1aμP2χ21+χ1.Substituting equations (46), (52), (61), and (67) into equation (39), and partially differentiating with respect to P, the partial derivative of g3 with respect to P is finally derived as follows:
(73)∂g3∂P=∂2J∂P∂fM+∂2J∂P2∂g2∂P−∂J∂P∂2g2∂P2⋅aM3/μ∂g1/∂fM−∂g2/∂fM∂g2/∂P2−∂J/∂P∂g2/∂g2⋅∂2g2∂P∂fM.Let us focus on the last term of the partial derivative ofg3. Partially differentiating equations (54), (58) ,and (70) with respect to fM, one can obtain that
(74)∂2r2∂fM2=2e2−e2cos2fM+ecosfM1+ecosfM2r2,∂2G∂fM2=sr12h2μP∂2r2∂fM2−2r1ph2r1h2sr12h22∂r2∂fM−r1ph22−r1h22r12+r1h24sr12h23r2,∂2F∂fM2=1r1P2∂r2∂fMr1ph2−r2r1h2−∂2r2∂fM2r1−r1h2.Partially differentiating equation (65) with respect to fM twice and introducing another set of notations are expressed as follows:
(75)JfMx3=∂2r2∂fM2h2x+2∂r2∂fM∂h2x∂fM−r2h2x1G−r2h2xG2∂2G∂fM2−2JfMx1G∂G∂fM−r1xJfMFG,JfMy3=∂2r2∂fM2⋅h2y+2∂r2∂fM∂h2y∂fM−r2h2y1G−r2h2yG2∂2G∂fM2−2JfMy1G∂G∂fM−r1yJfMFG,JfMz3=∂2r2∂fM2⋅h2z+2∂r2∂fM∂h2z∂fM−r2h2z1G−r2h2zG2∂2G∂fM2−2JfMz1G∂G∂fM−r1zJfMFG,where
(76)JfMFG=∂2F∂fM21G−FG2∂2G∂fM2−2G2∂F∂fM∂G∂fM+2FG3∂G∂fM2.Then, the second-order partial derivative of the cost functionJ with respect to fM is given by
(77)∂2J∂fM2=2JfMx22+JPx1JfMx3+JfMy22+JPy1JfMy3+JfMz22+JPz1JfMz3.Partially differentiating equation (60) with respect to fM, then one has
(78)∂2a∂fM2=P∂χ5∂fM1χ0−χ64∂χ51∂fM−2χ5χ51χ02+2χ64χ512χ0,where
(79)χ9=∂2r2∂fM2χ6−2∂r2∂fMr1ph2+r2r1h2,∂χ5∂fM=∂r2∂fM2+r2∂2r2∂fM22sr12h22−8r2r1h2r1ph2∂r2∂fM+2r22r1h22−r1ph22,∂χ51∂fM=r1−r2∂2r2∂fM2−χ9−∂r2∂fM22P2−2χ632−2r2χ6χ9+r1+2r2χ9−r2r1h2+2∂r2∂fM2χ62P.Partially differentiating equation (61), one can obtain the second-order partial derivative of g2 with respect to fM
(80)∂2g2∂fM2=∂2G∂fM2+32aμ∂2a∂fM2+12a∂a∂fM2arccosχ1−χ65Pχ2aμ∂a∂fM+aμχ91+1μP∂2a∂fM2χ7χ8+a∂2χ7∂fM2χ8+aχ7∂2χ8∂fM2+2∂a∂fM∂χ7∂fMχ8+2∂a∂fMχ7∂χ8∂fM+2a∂χ7∂fM∂χ8∂fM,where
(81)∂2χ7∂fM2=−r1h2r1ph2+r1ph2sr12h2r1h2sr12h2−χ7∂χ7∂fM,∂2χ8∂fM2=χ9P−∂2r2∂fM2.Furthermore, partially differentiating equation (46), the second-order partial derivative of g1 is expressed as follows:
(82)∂2g1∂fM2=eM1−eM2sinfM1+eMcosfM+eMsinEM.Then, the partial derivative ofg3 with respect to fM can be expressed as follows:
(83)∂g3∂fM=∂2J∂fM2+∂2J∂P∂fM∂g2∂P−∂J∂P∂2g2∂P∂fMaM3/μ∂g1/∂fM−∂g2/∂fM∂g2/∂P2+∂J/∂P∂g2/∂PaM3μ∂2g1∂fM2−∂2g2∂fM2.Now, the partial derivatives ofg3 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 variablesX=t,P,fMT for the MVIP. This method is based on a discovery that semilatus rectum, P, is insensitive to the variation of TOF in the computation process. In fact, the semilatus rectum P 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 P 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,Pis 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, t, and the true anomaly, fM. It should be noted that, according to equation (7), t can be analytically expressed as a function of fM. Therefore, the dimension of the search space can be further reduced to formulate a single-variable rooting problem. Let ΨfM denote the equality constraint on fM, which is expressed as follows:
(84)ΨfM=g2−aM3μg1.To determine the initial guess offM, a single-variable Newton-Raphson iteration equation is formulated as follows:
(85)fML+1=fML−ΨfMLΨ̇fML,where
(86)Ψ̇fM=dg2dfM−aM3μdg1dfM.It should be pointed out thatΨ̇fM can be calculated analytically from the derivation in Section 3.2.The iteration is finished if the difference betweentL and tL+1 is small enough, and the convergence criterion of the initialization method is given by
(87)tL+1−tL<εt.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εX, 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:
(88)XK+1−XK<εX.Otherwise, the iterative variablesX=t,P,fMT 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 rectumP reaches as high as 106, whereas the true anomaly fM 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 u and the corresponding convergence criterion is formulated as follows:
(89)uXK+1−uXK<εx,where u is the weight matrix.The procedure of implementing the proposed algorithm is included in the flowchart in Figure2. 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.Figure 2
Flowchart of the proposed method.


## 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 Table1 in detail.Table 1
Initial conditions for interception.
Group 1
Group 2
Group 3
rM0
5.9,0.84,2.8e+06
–7.4,2.8,–5.9e+06
5.2,3.7,–2.6e+06
vM0
–0.57,5.2,4.7e+03
–3.1,6.4,1.1e+03
0.3,6.2,2.9e+03
rE0
4.0,4.8,2.3e+06
0.22,6.2,−2.5e+06
1.7,6.3,–1.7e+06
vE0
Case 1
3.0,–3.0,5.0e+03
–5.0,0.0,–6.0e+03
6.0,2.0,1.0e+03
Case 2
1.0,2.0,5.0e+03
–4.0,3.0,–2.0e+03
3.0,4.0,3.0e+03
Case 3
–1.0,5.0,4.0e+03
–3.0,5.0,2.0e+03
1.0,5.0,3.0e+03


### 4.1. Convergent Characteristics

The three groups of simulation cases listed in Table1 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.Figure 3
The engagement geometry (a) and the iterative process (b) of Group 1.
(a)
(b)Figure 4
The engagement geometry and the iterative process of Group 2.
(a)
(b)Figure 5
The engagement geometry and the iterative process of Group 3.
(a)
(b)Table 2
Convergent property of the proposed method.
Convergent results of Group 1
Initialization method
Main iteration
L
fM
tL+1−tL
K
t,s
P,km
fM
uXK+1−uXK
Case 1
0
2.5000
/
0
479.99
4574.4
2.8318
/
JX=5.73e+05
1
2.8262
3.71e+02
1
448.06
5274.1
2.8056
7.01e–01
Δv=7.57e+02
2
2.8318
6.79e+00
2
453.18
5227.7
2.8097
4.69e–02
g1X=4.44e–16
3
2.8318
2.10e−02
3
453.25
5226.1
2.8098
1.54e–03
g2X=6.59e–12
4
2.8318
2.04e–07
4
453.26
5226.1
2.8098
3.97e–05
5
453.26
5226.1
2.8098
1.02e–06
6
453.26
5226.1
2.8098
2.61e–08
Case 2
0
3.2000
/
0
929.85
1755.9
3.1872
/
JX=1.04e+05
1
3.1818
2.33e+01
1
853.69
1777.5
3.1279
9.89e–02
Δv=3.22e+02
2
3.1868
6.39e+00
2
866.05
1793.6
3.1376
2.25e–02
g1X=4.44e–16
3
3.1872
5.71e–01
3
866.25
1794.0
3.1377
5.18e–04
g2X=7.52e–10
4
3.1872
4.61e–03
4
866.25
1794.0
3.1377
1.08e–06
5
866.25
1794.0
3.1377
1.75e–09
Case 3
0
4.0000
/
0
1664.93
2531.9
3.8067
/
JX=3.05e+05
1
3.8488
1.45e+02
1
1661.48
2526.6
3.8033
7.14e–03
Δv=5.52e+02
2
3.8090
4.05e+01
2
1661.36
2526.4
3.8032
2.31e–04
g1X=−7.22e–12
3
3.8067
2.41e+00
3
1661.36
2526.4
3.8032
8.38e–06
g2X=−7.66e–09
4
3.8067
7.90e–03
4
1661.36
2526.4
3.8032
3.04e–07As 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 Table2, two parts successfully converge to the required tolerance within six cycles. The maximal error for equality constraints, g1andg2, is less than 1e–8. The error norm of iterative parameters monotonously decreases and its maximal error is 3.04e–07. Therefore, the proposed method works well for all cases and has a high convergence rate.In addition, let us focus on the convergent process ofP 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 P 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
(90)εt=1e−02,εX=1e−06.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 Table1. The methods used for comparison include Zarchan’s method and the P-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 P-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.Table 3
Comparisons of the number of iteration.
Number of iteration (NOI)
Group 1
Group 2
Group 3
Case 1
Case 2
Case 3
Case 1
Case 2
Case 3
Case 1
Case 2
Case 3
Zarchan’s method
1676
1596
1737
1616
1564
1662
1614
1554
1657
P-iteration
174
90
118
144
56
88
212
160
126
Active set
18
30
14
20
17
17
19
15
19
Proposed method
6
5
4
4
3
6
5
6
6Table 4
Comparisons of the computation time.
Computation time (ms)
Group 1
Group 2
Group 3
Case 1
Case 2
Case 3
Case 1
Case 2
Case 3
Case 1
Case 2
Case 3
Zarchan’s method
4.46
7.86
7.03
8.14
6.87
4.18
10.04
10.50
9.45
P-iteration
2.78
4.01
4.07
3.37
5.83
1.86
4.13
6.81
7.20
Active set
32.42
26.86
12.50
22.60
16.59
17.48
18.10
12.37
15.86
Proposed method
0.090
0.078
0.067
0.090
0.103
0.116
0.078
0.089
0.090It can be seen from Table3 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 P-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 Table4. 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 P-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 Figure6. It should be noted that Zarchan’s method and the P-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.Figure 6
Efficiency comparison of four methods for solving the MVIP.


### 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 Table1. 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 Table5. 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.Table 5
Comparisons of simulation results.
Δv(m/s)
Group 1
Group 2
Group 3
Case 1
Case 2
Case 3
Case 1
Case 2
Case 3
Case 1
Case 2
Case 3
PN
1046
996
3575
802
721
2751
557
1244
1398
APN
1030
678
2673
753
371
2089
545
1097
1027
OTG
1025
592
2264
738
284
1738
543
1058
906
Proposed method
757
322
552
544
110
421
412
790
334
Minimum cost reduction (%)
26.1
45.6
75.6
26.3
61.3
75.8
24.1
25.3
63.1
Maximum cost reduction (%)
27.6
67.7
84.6
32.2
84.7
84.7
26.0
36.5
76.1


## 4.1. Convergent Characteristics

The three groups of simulation cases listed in Table1 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.Figure 3
The engagement geometry (a) and the iterative process (b) of Group 1.
(a)
(b)Figure 4
The engagement geometry and the iterative process of Group 2.
(a)
(b)Figure 5
The engagement geometry and the iterative process of Group 3.
(a)
(b)Table 2
Convergent property of the proposed method.
Convergent results of Group 1
Initialization method
Main iteration
L
fM
tL+1−tL
K
t,s
P,km
fM
uXK+1−uXK
Case 1
0
2.5000
/
0
479.99
4574.4
2.8318
/
JX=5.73e+05
1
2.8262
3.71e+02
1
448.06
5274.1
2.8056
7.01e–01
Δv=7.57e+02
2
2.8318
6.79e+00
2
453.18
5227.7
2.8097
4.69e–02
g1X=4.44e–16
3
2.8318
2.10e−02
3
453.25
5226.1
2.8098
1.54e–03
g2X=6.59e–12
4
2.8318
2.04e–07
4
453.26
5226.1
2.8098
3.97e–05
5
453.26
5226.1
2.8098
1.02e–06
6
453.26
5226.1
2.8098
2.61e–08
Case 2
0
3.2000
/
0
929.85
1755.9
3.1872
/
JX=1.04e+05
1
3.1818
2.33e+01
1
853.69
1777.5
3.1279
9.89e–02
Δv=3.22e+02
2
3.1868
6.39e+00
2
866.05
1793.6
3.1376
2.25e–02
g1X=4.44e–16
3
3.1872
5.71e–01
3
866.25
1794.0
3.1377
5.18e–04
g2X=7.52e–10
4
3.1872
4.61e–03
4
866.25
1794.0
3.1377
1.08e–06
5
866.25
1794.0
3.1377
1.75e–09
Case 3
0
4.0000
/
0
1664.93
2531.9
3.8067
/
JX=3.05e+05
1
3.8488
1.45e+02
1
1661.48
2526.6
3.8033
7.14e–03
Δv=5.52e+02
2
3.8090
4.05e+01
2
1661.36
2526.4
3.8032
2.31e–04
g1X=−7.22e–12
3
3.8067
2.41e+00
3
1661.36
2526.4
3.8032
8.38e–06
g2X=−7.66e–09
4
3.8067
7.90e–03
4
1661.36
2526.4
3.8032
3.04e–07As 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 Table2, two parts successfully converge to the required tolerance within six cycles. The maximal error for equality constraints, g1andg2, is less than 1e–8. The error norm of iterative parameters monotonously decreases and its maximal error is 3.04e–07. Therefore, the proposed method works well for all cases and has a high convergence rate.In addition, let us focus on the convergent process ofP 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 P 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
(90)εt=1e−02,εX=1e−06.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 Table1. The methods used for comparison include Zarchan’s method and the P-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 P-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.Table 3
Comparisons of the number of iteration.
Number of iteration (NOI)
Group 1
Group 2
Group 3
Case 1
Case 2
Case 3
Case 1
Case 2
Case 3
Case 1
Case 2
Case 3
Zarchan’s method
1676
1596
1737
1616
1564
1662
1614
1554
1657
P-iteration
174
90
118
144
56
88
212
160
126
Active set
18
30
14
20
17
17
19
15
19
Proposed method
6
5
4
4
3
6
5
6
6Table 4
Comparisons of the computation time.
Computation time (ms)
Group 1
Group 2
Group 3
Case 1
Case 2
Case 3
Case 1
Case 2
Case 3
Case 1
Case 2
Case 3
Zarchan’s method
4.46
7.86
7.03
8.14
6.87
4.18
10.04
10.50
9.45
P-iteration
2.78
4.01
4.07
3.37
5.83
1.86
4.13
6.81
7.20
Active set
32.42
26.86
12.50
22.60
16.59
17.48
18.10
12.37
15.86
Proposed method
0.090
0.078
0.067
0.090
0.103
0.116
0.078
0.089
0.090It can be seen from Table3 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 P-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 Table4. 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 P-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 Figure6. It should be noted that Zarchan’s method and the P-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.Figure 6
Efficiency comparison of four methods for solving the MVIP.


## 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 Table1. 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 Table5. 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.Table 5
Comparisons of simulation results.
Δv(m/s)
Group 1
Group 2
Group 3
Case 1
Case 2
Case 3
Case 1
Case 2
Case 3
Case 1
Case 2
Case 3
PN
1046
996
3575
802
721
2751
557
1244
1398
APN
1030
678
2673
753
371
2089
545
1097
1027
OTG
1025
592
2264
738
284
1738
543
1058
906
Proposed method
757
322
552
544
110
421
412
790
334
Minimum cost reduction (%)
26.1
45.6
75.6
26.3
61.3
75.8
24.1
25.3
63.1
Maximum cost reduction (%)
27.6
67.7
84.6
32.2
84.7
84.7
26.0
36.5
76.1


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



---

*Source: 8502870-2019-09-16.xml*