Damped Iterative Explicit Guidance for Multistage Rockets with Thrust Drop Faults

Abstract

A damped iterative explicit guidance (DIEG) algorithm is proposed to address the problem of the insufficient convergence of classical explicit guidance methods in the event of thrust drop faults in multistage rockets. Based on the iterative guidance mode (IGM) and powered explicit guidance (PEG), this method is enhanced in three aspects: (1) an accurate transversality condition is derived and applied in the dimension-reduction framework instead of using a simplified assumption; (2) the Gauss–Legendre quadrature formula (GLQF) is adopted to increase the accuracy of the method by addressing the issue of excessive errors in calculating thrust integration using linearization methods based on a small quantity assumption under fault conditions; and (3) a damping factor for solving the time-to-go is introduced to avoid the chattering phenomenon and enhance convergence. A numerical simulation was conducted in single- and multistage mission scenarios by gradually reducing the engine thrust to compare the performance of DIEG and PEG. The results show that DIEG has a much larger convergence range than PEG and has fuel optimality similar to that of the optimization method in most fault scenarios. Finally, the robustness of DIEG under various deviations is verified through Monte Carlo simulation.

1. Introduction

Improving the reliability of launch vehicles, especially enhancing their fault tolerance, is an important direction for the development of launch vehicle technology. Based on their severity, the faults can be classified into small-scale faults with the ability to enter the original target orbit, medium-scale faults with the ability to enter the downgraded orbit, and catastrophic faults with the inability to enter any orbit [1]. The onboard guidance algorithm should enable the rocket to enter its original or downgrade orbit with sufficient remaining fuel. Therefore, the guidance algorithm must be characterized by strong optimality and wide adaptability to faults.

Early rockets used perturbation guidance that tracked standard trajectories, which were simple to implement, but simultaneously satisfying multiple constraints was difficult, and resistance to large deviation was weak, so this method was unsuitable for fault situations. Since the 1960s, explicit guidance algorithms [2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22] that do not rely on standard trajectories have increasingly been applied for the extra-atmospheric guidance of launch vehicles. These algorithms derive the optimal explicit expression of the optimal control variable with respect to time; under various simplified assumptions, the command angle is iteratively calculated by reducing the dimensions based on the current flight state and target orbit constraints. Therefore, these algorithms provide advantages in accuracy, optimality, and robustness compared with the perturbation guidance method. Such algorithms can be divided into two categories: iterative guidance mode (IGM), originating from the Saturn V rocket [2,3,4,5], and powered explicit guidance (PEG), originating from space shuttles [6,7,8,9,10]. IGM simplifies the command angle into the sum of two parts, namely, the primary part in a constant form determined by velocity constraints and the secondary part in a constant plus linear form determined by position constraints. The time-to-go (t_go) and the main part are iteratively solved based on velocity constraints first; second, the complete command angle that satisfies both the velocity and position constraints is calculated by assuming that the quantity of the secondary part is small. In a single guidance period, IGM solves the command angle through an open-loop method without predicting or correcting the entry point. Therefore, the solution may contain considerable errors in the early stage of guidance relative to the optimal solution, which limits its application range. PEG uses a three-dimensional thrust vector as the control variable, and the orthogonal assumption is adopted to obtain the linear tangent form of the command angle. The velocity to be gained (v_go) is used as the iterative variable to predict and correct the injection point’s state, which improves thrust integration and gravity integration methods. This algorithm provides advantages in terms of optimality and adaptability compared with IGM owing to these features.

IGM and PEG have been modified and enhanced to suit a wide range of mission needs since their introduction [11,12,13,14,15,16,17]. Researchers [11,12,13] have proposed higher-precision thrust integration methods and adopted linear differential correction to accelerate convergence. Song [14] improved the initial estimate of v_go in PEG, thereby enhancing the convergence of the original algorithm when approaching orbit. Researchers [15,16] have enhanced multiple aspects of the original algorithm, including limiting control variables and scaling identity Jacobian, in order to apply PEG to the Space Launch System (SLS), to improve its performance in long burn arcs. Inspired by PEG, IGM modifications have mainly focused on introducing higher-order terms in thrust integration as well as the prediction and correction of the position of the injection point [14,19,20,21,22]. The optimality of enhanced IGM is similar to that of IGM [14] and does not require the iteration of v_go, thereby reducing the risk of nonconvergence. However, the above enhancement algorithms retain the simplified transversality condition and the linearized thrust integration method based on the small quantity assumption. During a thrust drop fault, the simplified assumptions above may no longer hold due to the increase in position correction, requiring increased fuel consumption and potentially causing the algorithm to fail to converge. Therefore, the existing explicit guidance algorithms must be further improved to enhance their fault tolerance capability.

Optimal guidance (OPGUID) [23,24,25,26,27,28] is another type of guidance method that is based on optimal control theory and the assumption of a linear gravitational field to obtain analytical expressions for the control variables and explicit integrals of the state variables with respect to time. The shooting method is used to solve the initial values of the costate variable. This method does not introduce any simplifying assumptions other than gravity; therefore, the optimality is stronger than that of explicit guidance algorithms. However, practical application is challenging due to the sensitivity of this method to the initial estimate of the costate variable [29]. Scholars are increasingly using optimization algorithms with the improvement in the capability of onboard computers [30,31,32,33,34], which are represented by convex optimization for rocket guidance. These algorithms are highly accurate and handling constraints is easy, so they are suitable for various launch missions. However, various nonconvex terms, especially the constraint of nonconvex terminal orbital elements, need sequence linearization to achieve convexity. Guaranteeing the convergence of this process is difficult [35], which restricts practical applications in ascending guidance missions. In addition, some scholars are exploring the use of machine learning methods [36,37,38,39,40] for rocket guidance, but the noninterpretability of deep neural networks can create unpredictable risks for rockets. In summary, explicit guidance algorithms will remain the preferred choice for engineering applications for the foreseeable future.

On the basis of the above studies, we developed a damped iterative explicit guidance (DIEG) algorithm to enhance the adaptability of multistage rockets in the thrust drop fault situation. The main contributions of this study are as follows:

(1) We derived the accurate transversality condition of the optimal solution with terminal orbital element constraints and applied them to the dimension-reduction iterative framework of explicit guidance algorithms.

(2) We applied the Gauss–Legendre quadrature formula (GLQF) to the multistage thrust integration calculation, thereby avoiding the problems caused by the small quantity assumption’s failure in the linearization methods used in previous explicit guidance algorithms.

(3) A damping factor was introduced to the iterative calculation of the time-to-go to avoid chattering and substantially improve convergence.

The remaining sections of this article are organized as follows: In Section 2, the optimal control problem for a multistage rocket entering orbit is described, and the optimality conditions including transversality conditions are derived. Section 3 outlines the method of using the GLQF to calculate the multistage rocket thrust integral. Section 4 outlines the iterative process of the DIEG algorithm, with Section 4.1 focusing on the method of introducing the damping factor to solve the time-to-go. In Section 5, the numerical simulations conducted to compare and verify the effectiveness of the various improvements in the DIEG algorithm, as well as its robustness, are described.

2. Optimal Control Problem for Multistage Rocket Entering Orbit

2.1. Optimal Control Problem Model

The characteristics of the studied multistage guidance problem are as follows: (1) The rocket undergoes multiple burn phases from the current moment to the injection point, and the thrust magnitude of each burn phase may be different but cannot be throttled. (2) A coast phase may occur between any two adjacent burn phases, and the coast time is preset. (3) The fuel consumed in each phase, except for the final phase, is preset, so the minimum fuel consumption in the final phase is the indicator for evaluating the optimality of the guidance algorithm.

First, define the reference injection point coordinate system O-XYZ as shown in Figure 1. The coordinate origin is taken as the center of the Earth; the OY axis points to the reference injection point and is fixed in inertial space; the OX axis is in the orbital plane, perpendicular to the OY axis and pointing in the direction of satellite movement; and the OZ, OX, and OY axes form a right-handed coordinate system. In Figure 1, O-X_IY_IZ_I is the geocentric equatorial inertial coordinate system; i* is the inclination angle of the target orbit; w* is the argument of perigee; Ω* is the ascending node right ascension; and f* is the true anomaly of the reference injection point.

Establish the optimal control problem model for the multistage rocket entering orbit in O-XYZ. The equation of the state is

$$\left\{\begin{array}{l}{\dot{\mathit{r}}}_i={\mathit{v}}_i\\ {\dot{\mathit{v}}}_i=\mathit{g}+{\displaystyle \frac{T_i}{m}}\cdot {\mathit{u}}_i\\ {\dot{m}}_i=-m_{si}, i=1,2,\cdots ,n\end{array}\right.$$

(1)

where subscript i represents the i-th flight phase; r represents the rocket position vector; v represents the velocity vector; m represents mass; g represents the gravitational vector; u represents the unit thrust vector; T represents the constant thrust magnitude; and m_s represents the engine flow rate per second. When i is the coasting phase, T_i = 0, and m_si = 0. Adopt the assumption of a linear gravitational field,

$$\mathit{g}=-{\omega }^2\mathit{r}$$

(2)

where $\omega$ is the average angular rate, taken as the angular rate of the circular orbit at the average height of current position r₀ and the injection point position r_f; that is,

$$\begin{array}{l}\overline{r}={\displaystyle \frac{‖{\mathit{r}}_0‖+‖{\mathit{r}}_{\mathrm{f}}‖}{2}}\\ \omega =\sqrt{{\displaystyle \frac{\mu }{{\overline{r}}^3}}}\end{array}$$

(3)

where $\mu$ is the gravitational constant of the Earth.

The minimum fuel consumption in the final phase is equivalent to the shortest flight time due to the constant duration of each phase (including the coasting phase), except for the final phase. Therefore, the cost function of this optimal control problem is

$$\min {\displaystyle {\int }_0^{t_{\mathrm{go}}}1}\mathrm{d}t$$

(4)

where $t_{\mathrm{go}}$ is the time-to-go, which is the total flight time from the current moment until the moment of injection. The constraint on the control variables is

$$‖{\mathit{u}}_i‖=1$$

(5)

The terminal constraints include five orbital elements, namely, the semimajor axis a, eccentricity e, orbital inclination i, ascending node right ascension Ω, and perigee argument w, which can be converted into

$$\begin{array}{l}{C}_1=v_{\mathrm{f}}^2-{\displaystyle \frac{2\mu }{r_{\mathrm{f}}}}+{\displaystyle \frac{\mu }{a^{\ast }}}=0\\ C_2=\left({\mathit{r}}_{\mathrm{f}}\times {\mathit{v}}_{\mathrm{f}}\right)\cdot \left({\mathit{r}}_{\mathrm{f}}\times {\mathit{v}}_{\mathrm{f}}\right)-\mu \cdot a^{\ast }\left(1-e^{\ast 2}\right)=0\\ C_3=\left[\left(v^2-{\displaystyle \frac{\mu }{r_{\mathrm{f}}}}\right){\mathit{r}}_{\mathrm{f}}-\left({\mathit{r}}_{\mathrm{f}}\cdot {\mathit{v}}_{\mathrm{f}}\right){\mathit{v}}_{\mathrm{f}}\right]\cdot {\displaystyle \frac{\mathit{N}}{\mu }}-\cos \left(w^{\ast }\right)e^{\ast }=0\\ C_4=v_{z\mathrm{f}}=0\\ C_5=r_{z\mathrm{f}}=0\end{array}$$

(6)

where subscript f represents the terminal value; subscripts x, y, and z represent the components of a vector along the three axes of O-XYZ, respectively; a* is the semimajor axis of the target orbit; e* is the eccentricity of the target orbit; and $\mathit{N}$ represents the unit vector of the intersection line (pointing from the center of the earth to the ascending node) between the target orbital and the equatorial planes.

2.2. Solving Optimal Control Problem

The Hamiltonian function of the optimal control problem in Section 2.1 is

$$H=-1+{\mathit{\lambda }}_r\cdot \mathit{v}+{\mathit{\lambda }}_{\mathit{v}}\cdot \left(-{\mathit{\omega }}^2\mathit{r}+{\displaystyle \frac{T}{m}}\cdot \mathit{u}\right)+{\mathit{\lambda }}_m\cdot \left(-m_{\mathrm{s}}\right)$$

(7)

where λ_r, λ_v, and λ_m are costate variables. We did not study the switching function such as OPGUID because the segmentation time was preset. According to the Pontryagin maximum principle, the control variable is solved as

$$\mathit{u}={\displaystyle \frac{{\mathit{\lambda }}_v}{‖{\mathit{\lambda }}_v‖}}$$

(8)

and the costate variable are [8]

$$\left\{\begin{array}{l}{\mathit{\lambda }}_v={\mathit{\lambda }}_{v0}\cos \left(\omega t\right)-{\displaystyle \frac{{\mathit{\lambda }}_{r0}}{\omega }}\sin \left(\omega t\right)\\ {\mathit{\lambda }}_r=\omega {\mathit{\lambda }}_{v0}\sin \left(\omega t\right)+{\mathit{\lambda }}_{r0}\cos \left(\omega t\right)\end{array}\right.$$

(9)

Let $\mathit{\lambda }={\mathit{\lambda }}_v$ and $\dot{\mathit{\lambda }}=-{\mathit{\lambda }}_r$; by introducing the control parameter K, where ${\mathit{\lambda }}_K, {\dot{\mathit{\lambda }}}_K$, respectively, represent the values of $\mathit{\lambda }, \dot{\mathit{\lambda }}$ at time K, Equation (9) can be transformed into

$$\left\{\begin{array}{l}\mathit{\lambda }={\mathit{\lambda }}_K\cos \left[\omega \left(t-K\right)\right]+{\displaystyle \frac{{\dot{\mathit{\lambda }}}_K}{\omega }}\sin \left[\omega \left(t-K\right)\right]\\ \dot{\mathit{\lambda }}=-\omega {\mathit{\lambda }}_K\sin \left[\omega \left(t-K\right)\right]+{\dot{\mathit{\lambda }}}_K\cos \left[\omega \left(t-K\right)\right]\end{array}\right.$$

(10)

In addition, the expression of the transversality conditions of the solution to the optimal control problem is

$$\left\{\begin{array}{l}{\mathit{\lambda }}_{\mathrm{f}}=-{\displaystyle \sum _{i=1}^5{\eta }_i\frac{\partial C_i}{\partial {\mathit{v}}_{\mathrm{f}}}}\\ {\dot{\mathit{\lambda }}}_{\mathrm{f}}={\displaystyle \sum _{i=1}^5{\eta }_i\frac{\partial C_i}{\partial {\mathit{r}}_{\mathrm{f}}}}\end{array}\right.$$

(11)

and

$$H_{\mathrm{f}}=0$$

(12)

where ${\eta }_i$ represents the Lagrange multiplier. The constraint in Equation (12) can be replaced by Equation (13) [24]:

$$‖{\mathit{\lambda }}_K‖=1$$

(13)

By substituting Equation (6) into Equation (11), we obtain

$${\mathit{\lambda }}_{\mathrm{f}}=-2{\eta }_1{\mathit{v}}_{\mathrm{f}}-{\eta }_2\left({\mathit{r}}_{\mathrm{f}}\times {\mathit{v}}_{\mathrm{f}}\right)\times {\mathit{r}}_{\mathrm{f}}-{\eta }_4{e}_3-{\eta }_3\left[2{\mathit{v}}_{\mathrm{f}}\left({\mathit{r}}_{\mathrm{f}}\cdot \mathrm{N}\right)-{\mathit{r}}_{\mathrm{f}}\left({\mathit{v}}_{\mathrm{f}}\cdot \mathit{N}\right)-\left({\mathit{r}}_{\mathrm{f}}\cdot {\mathit{v}}_{\mathrm{f}}\right)\mathit{N}\right]$$

(14)

$${\dot{\mathit{\lambda }}}_{\mathrm{f}}=2{\eta }_1\mu {\displaystyle \frac{{\mathit{r}}_{\mathrm{f}}}{r_{\mathrm{f}}^3}}+{\eta }_2{\mathit{v}}_{\mathrm{f}}\times \left({\mathit{r}}_{\mathrm{f}}\times {\mathit{v}}_{\mathrm{f}}\right)+{\eta }_5{e}_3+{\eta }_3\left[\left(v_{\mathrm{f}}^2-{\displaystyle \frac{\mu }{r_{\mathrm{f}}}}\right)\mathit{N}+\mu {\displaystyle \frac{{\mathit{r}}_{\mathrm{f}}}{r_{\mathrm{f}}^3}}\left({\mathit{r}}_{\mathrm{f}}\cdot \mathit{N}\right)-\left({\mathit{v}}_{\mathrm{f}}\cdot \mathit{N}\right){\mathit{v}}_{\mathrm{f}}\right]$$

(15)

where $e_3={\left[0,0,1\right]}^{\mathrm{T}}$. Simultaneously multiply both sides of Equation(14) by $-\mu {\displaystyle \frac{{\mathit{r}}_{\mathrm{f}}}{r_{\mathrm{f}}^3}}$ to obtain

$$-\mu {\displaystyle \frac{{\mathit{r}}_{\mathrm{f}}}{r_{\mathrm{f}}^3}}\cdot {\mathit{\lambda }}_{\mathrm{f}}=2{\eta }_1\mu {\displaystyle \frac{{\mathit{r}}_{\mathrm{f}}\cdot {\mathit{v}}_{\mathrm{f}}}{r_{\mathrm{f}}^3}}+{\eta }_3\mu {\displaystyle \frac{1}{r_{\mathrm{f}}^3}}\left[\left({\mathit{v}}_{\mathrm{f}}\cdot {\mathit{r}}_{\mathrm{f}}\right)\left({\mathit{r}}_{\mathrm{f}}\cdot \mathit{N}\right)-r_{\mathrm{f}}^2\left({\mathit{v}}_{\mathrm{f}}\cdot \mathit{N}\right)\right]$$

(16)

Simultaneously multiply both sides of Equation (15) by ${\mathit{v}}_{\mathrm{f}}$ to obtain

$${\dot{\mathit{\lambda }}}_{\mathrm{f}}\cdot {\mathit{v}}_{\mathrm{f}}=2{\eta }_1\mu {\displaystyle \frac{{\mathit{r}}_{\mathrm{f}}\cdot {\mathit{v}}_{\mathrm{f}}}{r_{\mathrm{f}}^3}}+{\eta }_3\mu {\displaystyle \frac{1}{r_{\mathrm{f}}^3}}\left[-r_{\mathrm{f}}^2\left(\mathit{N}\cdot {\mathit{v}}_{\mathrm{f}}\right)+\left({\mathit{r}}_{\mathrm{f}}\cdot \mathit{N}\right)\left({\mathit{v}}_{\mathrm{f}}\cdot {\mathit{r}}_{\mathrm{f}}\right)\right]$$

(17)

The Lagrange multiplier can be eliminated by subtracting Equation (17) from Equation (16), and the final expression for the transversality condition is obtained as follows:

$$-\mu {\displaystyle \frac{{\mathit{r}}_{\mathrm{f}}}{r_{\mathrm{f}}^3}}\cdot {\mathit{\lambda }}_{\mathrm{f}}={\dot{\mathit{\lambda }}}_{\mathrm{f}}\cdot {\mathit{v}}_{\mathrm{f}}$$

(18)

Explicit guidance methods usually use simplified assumptions instead of accurate transversality conditions (18). For example, PEG uses orthogonality assumptions [7]:

$${\mathit{\lambda }}_K\cdot {\dot{\mathit{\lambda }}}_K=0$$

(19)

whereas IGM uses three velocity and two position component constraints instead of five orbital element constraints [19], resulting in

$${\dot{\lambda }}_{\mathrm{f}X}=0$$

(20)

Our proposed method thus retains an accurate transversality condition and theoretically increases the algorithm’s optimality.

3. High-Precision Thrust Integration Using GLQF

The explicit guidance algorithm requires the use of thrust integration for solving the control parameters and predicting the terminal states, so the accuracy of thrust integration directly impacts the accuracy and optimality of the explicit guidance algorithm. Both IGM and PEG use linearization methods to integrate thrust, but a large position deviation increases $‖{\dot{\mathit{\lambda }}}_K‖$ when a thrust fault occurs. Therefore, the assumption used in the linearization process that the modulus of ${\dot{\mathit{\lambda }}}_K\left(t_{\mathrm{go}}-K\right)$ is small may not be met, leading to substantial thrust integration errors. In response to this issue, we adopted the GLQF method to integrate thrust.

First, the thrust integral is defined as follows:

$${\mathit{v}}_{\mathrm{thrust}}={\displaystyle {\int }_0^{t_{\mathrm{go}}}\frac{T}{m}\mathit{u}}\mathrm{d}t=L_{\mathrm{T}}{\mathit{\lambda }}_K+J_{\mathrm{T}}{\dot{\mathit{\lambda }}}_K$$

(21)

$${\mathit{r}}_{\mathrm{thrust}}={\displaystyle {\int }_0^{t_{\mathrm{go}}}{\displaystyle {\int }_0^s\frac{T}{m}\mathit{u}}}\mathrm{d}t\mathrm{d}s=S_{\mathrm{T}}{\mathit{\lambda }}_K+Q_{\mathrm{T}}{\dot{\mathit{\lambda }}}_K$$

(22)

where L_T, J_T, S_T, and Q_T are the full integration coefficients of thrust, expressed as

$$\begin{array}{l}{L}_{\mathrm{T}}={\displaystyle {\int }_0^{t_{\mathrm{go}}}\frac{\cos \left[\omega \left(t-K\right)\right]}{d_{\mathrm{en}}}\frac{T}{m}}\mathrm{d}t\\ J_{\mathrm{T}}={\displaystyle {\int }_0^{t_{\mathrm{go}}}\frac{\sin \left[\omega \left(t-K\right)\right]}{\omega d_{\mathrm{en}}}\frac{T}{m}}\mathrm{d}t\\ S_{\mathrm{T}}={\displaystyle {\int }_0^{t_{\mathrm{go}}}{\displaystyle {\int }_0^s\frac{\cos \left[\omega \left(t-K\right)\right]}{d_{\mathrm{en}}}\frac{T}{m}}}\mathrm{d}t\mathrm{d}s\\ Q_{\mathrm{T}}={\displaystyle {\int }_0^{t_{\mathrm{go}}}{\displaystyle {\int }_0^s\frac{\sin \left[\omega \left(t-K\right)\right]}{\omega d_{\mathrm{en}}}\frac{T}{m}}}\mathrm{d}t\mathrm{d}s\end{array}$$

(23)

where

$$\begin{array}{l}{d}_{\mathrm{en}}=\sqrt{1+\left({\displaystyle \frac{D_1-{\omega }^2}{{\omega }^2}}\right){\sin }^2\left[\omega \left(t-K\right)\right]+{\displaystyle \frac{D_2}{\omega }}\sin \left[2\omega \left(t-K\right)\right]}\\ D_1={‖{\dot{\mathit{\lambda }}}_K‖}^2\\ D_2={\mathit{\lambda }}_K\cdot {\dot{\mathit{\lambda }}}_K\end{array}$$

(24)

Before solving the control parameters, the values from the previous iteration are used for D₁ and D₂, and the algorithm is initialized with $D_1={\omega }^2$ and $D_2=0$.

The mass-to-mass flow rate ratio τ_i in phase i is defined as follows:

$${\tau }_i=m_{0i}/m_{\mathrm{s}i}$$

(25)

where m_0i is the initial mass in phase i. When i = 1, τ_i is calculated from the apparent acceleration’s modulus a_c0 at the current moment and the specific impulse I_sp1 at the current phase:

$${\tau }_1=I_{\mathrm{sp}1}/a_{\mathrm{c}0}$$

(26)

Let t_sum(i) be the total time from the current moment to the end of phase i (note that t_sum(0) = 0 and t_sum(n) = t_go). We have

$${\displaystyle \frac{T_i}{m_i}}={\displaystyle \frac{I_{\mathrm{sp}i}}{{\tau }_i-\left(t-t_{\mathrm{sum}\left(i-1\right)}\right)}}$$

(27)

where I_spi is the specific impulse in phase i. Let t_bi be the estimated burn time remaining in phase i, with $t_{\mathrm{b}n}=t_{\mathrm{go}}-t_{\mathrm{sum}\left(n-1\right)}$ in the final phase and preset values in all other phases. t_b1 needs to be updated when a fault occurs to $t_{\mathrm{b}1}/\left(1-\eta \right)$, where η is the thrust drop ratio. Then, the partial thrust integration coefficients L_i, J_i, S_i, and Q_i are defined in phase i. When i is the coast phase, $L_i=0$, $J_i=0$, $S_i=0$, and $Q_i=0$. When i is the burn phase,

$$L_i={\displaystyle {\int }_0^{t_{\mathrm{b}i}}\frac{\cos \left[\omega \left(t+t_{\mathrm{sum}\left(i-1\right)}-K\right)\right]}{d_{\mathrm{en}i}}\frac{I_{\mathrm{sp}i}}{{\tau }_i-t}}\mathrm{d}t={\displaystyle {\int }_0^{t_{\mathrm{b}i}}{f}_{Li}\left(t\right)}\mathrm{d}t$$

(28)

$$J_i={\displaystyle {\int }_0^{t_{\mathrm{b}i}}\frac{\sin \left[\omega \left(t+t_{\mathrm{sum}\left(i-1\right)}-K\right)\right]}{\omega d_{\mathrm{en}i}}\frac{I_{\mathrm{sp}i}}{{\tau }_i-t}}\mathrm{d}t={\displaystyle {\int }_0^{t_{\mathrm{b}i}}{f}_{Ji}\left(t\right)}\mathrm{d}t$$

(29)

$$S_i={\displaystyle {\int }_0^{t_{\mathrm{b}i}}{\displaystyle {\int }_0^s\frac{\cos \left[\omega \left(t+t_{\mathrm{sum}\left(i-1\right)}-K\right)\right]}{d_{\mathrm{en}i}}\frac{I_{\mathrm{sp}i}}{{\tau }_i-t}}}\mathrm{d}t\mathrm{d}s={\displaystyle {\int }_0^{t_{\mathrm{b}i}}{\displaystyle {\int }_0^s{f}_{Si}\left(t\right)}}\mathrm{d}t\mathrm{d}s$$

(30)

$$Q_i={\displaystyle {\int }_0^{t_{\mathrm{b}i}}{\displaystyle {\int }_0^s\frac{\sin \left[\omega \left(t+t_{\mathrm{sum}\left(i-1\right)}-K\right)\right]}{\omega d_{\mathrm{en}i}}\frac{I_{\mathrm{sp}i}}{{\tau }_i-t}}}\mathrm{d}t\mathrm{d}s={\displaystyle {\int }_0^{t_{\mathrm{b}i}}{\displaystyle {\int }_0^s{f}_{Qi}\left(t\right)}}\mathrm{d}t\mathrm{d}s$$

(31)

where

$$d_{\mathrm{en}i}=\sqrt{1+\left({\displaystyle \frac{D_1-{\omega }^2}{{\omega }^2}}\right){\sin }^2\left[\omega \left(t+t_{\mathrm{sum}\left(i-1\right)}-K\right)\right]+{\displaystyle \frac{D_2}{\omega }}\sin \left[2\omega \left(t+t_{\mathrm{sum}\left(i-1\right)}-K\right)\right]}$$

(32)

L_i, J_i, S_i, and Q_i can be quickly and accurately calculated using the GLQF method because they are all explicit integrals with respect to time, which is described in detail below.

For integrals L_i and J_i, first let

$$\alpha ={\displaystyle \frac{2t-t_{\mathrm{b}i}}{t_{\mathrm{b}i}}}$$

(33)

Then, we obtain

$$t={\displaystyle \frac{1}{2}}\left(\alpha +1\right)t_{\mathrm{b}i}$$

(34)

and

$$\begin{array}{l}{L}_i={\displaystyle {\int }_{-1}^1{g}_{Li}\left(\alpha \right)}\mathrm{d}\alpha \\ J_i={\displaystyle {\int }_{-1}^1{g}_{Ji}\left(\alpha \right)}\mathrm{d}\alpha \end{array}$$

(35)

where

$$\begin{array}{l}{g}_{Li}\left(\alpha \right)=f_{Li}\left({\displaystyle \frac{1}{2}}\left(\alpha +1\right)t_{\mathrm{b}i}\right)\cdot {\displaystyle \frac{1}{2}}{t}_{\mathrm{b}i}\\ g_{Ji}\left(\alpha \right)=f_{Ji}\left({\displaystyle \frac{1}{2}}\left(\alpha +1\right)t_{\mathrm{b}i}\right)\cdot {\displaystyle \frac{1}{2}}{t}_{\mathrm{b}i}\end{array}$$

(36)

According to GLQF, the value of integral (35) can be obtained as

$$\begin{array}{l}{L}_i\approx {\displaystyle \sum _{j=0}^n{w}_j^{\prime }{g}_{Li}\left(s_j\right)}\\ J_i\approx {\displaystyle \sum _{j=0}^n{w}_j^{\prime }{g}_{Ji}\left(s_j\right)}\end{array}$$

(37)

where s_j is a node of the Legendre polynomial, and $w_j^{\prime }$ is the corresponding weight.

For double integrals S_i and Q_i, let

$$\begin{array}{l}\alpha ={\displaystyle \frac{2s-t_{\mathrm{b}i}}{t_{\mathrm{b}i}}}\\ \beta ={\displaystyle \frac{2t-s}{s}}\end{array}$$

(38)

We thus obtain

$$\begin{array}{l}s={\displaystyle \frac{1}{2}}\left(\alpha +1\right)t_{\mathrm{b}i}\\ t={\displaystyle \frac{1}{4}}\left(\alpha +1\right)\left(\beta +1\right)t_{\mathrm{b}i}\end{array}$$

(39)

The Jacobian determinant is

$$\left|{\displaystyle \frac{\partial \left(t,s\right)}{\partial \left(\beta ,\alpha \right)}}\right|={\displaystyle \frac{1}{8}}\left(\alpha +1\right)t_{\mathrm{b}i}^2$$

(40)

As such, the integration of S_i and Q_i can be converted into

$$\begin{array}{l}{S}_i={\displaystyle {\int }_{-1}^1\mathrm{d}\alpha {\displaystyle {\int }_{-1}^1{g}_{Si}\left(\alpha ,\beta \right)}}\mathrm{d}\beta \\ Q_i={\displaystyle {\int }_{-1}^1\mathrm{d}\alpha {\displaystyle {\int }_{-1}^1{g}_{Qi}\left(\alpha ,\beta \right)}}\mathrm{d}\beta \end{array}$$

(41)

where

$$\begin{array}{l}{g}_{Si}\left(\alpha ,\beta \right)=f_{Si}\left({\displaystyle \frac{1}{4}}\left(\alpha +1\right)\left(\beta +1\right)t_{\mathrm{b}i}\right)\cdot {\displaystyle \frac{1}{8}}\left(\alpha +1\right)t_{\mathrm{b}i}^2\\ g_{Qi}\left(\alpha ,\beta \right)=f_{Qi}\left({\displaystyle \frac{1}{4}}\left(\alpha +1\right)\left(\beta +1\right)t_{\mathrm{b}i}\right)\cdot {\displaystyle \frac{1}{8}}\left(\alpha +1\right)t_{\mathrm{b}i}^2\end{array}$$

(42)

The value of integral (41) is

$$\begin{array}{l}{S}_i\approx {\displaystyle \sum _{j=0}^n{\displaystyle \sum _{k=0}^n{w}_j^{\prime }{w}_k^{\prime }{g}_{Si}\left(s_j,s_k\right)}}\\ Q_i\approx {\displaystyle \sum _{j=0}^n{\displaystyle \sum _{k=0}^n{w}_j^{\prime }{w}_k^{\prime }{g}_{Qi}\left(s_j,s_k\right)}}\end{array}$$

(43)

The accuracy of the GLQF method depends on the number of nodes. In this study, only six nodes were used to obtain high-precision integration values without consuming excessive computational resources.

Next, calculate the full integration coefficient of thrust, where

$$L_T={\displaystyle \sum _{i=1}^n{L}_i},J_T={\displaystyle \sum _{i=1}^n{J}_i}$$

(44)

S_i and Q_i need to be modified as follows before calculating the full integration coefficient:

$$\begin{array}{l}{S}_i=S_i+t_{bi}\cdot L_{Ti}\\ Q_i=Q_i+t_{bi}\cdot J_{Ti}\end{array}$$

(45)

where

$$L_{Ti}={\displaystyle \sum _{j=1}^i{L}_j},J_{Ti}={\displaystyle \sum _{j=1}^i{J}_j}$$

(46)

Finally, sum S_i and Q_i to obtain the full integration coefficient.

$$S_T={\displaystyle \sum _{i=1}^n{S}_i},Q_T={\displaystyle \sum _{i=1}^n{Q}_i}$$

(47)

4. Damped Iterative Explicit Guidance

The DIEG proposed in this article combines and enhances the iterative processes of IGM and PEG. The specific details of the algorithm are provided below.

4.1. Damped Iteration Method for Time-to-Go

First, define the velocity to be gained ${\mathit{v}}_{\mathrm{gn}}$ and position to be gained ${\mathit{r}}_{\mathrm{gn}}$:

$${\mathit{v}}_{\mathrm{gn}}={\mathit{v}}_{\mathrm{d}}-{\mathit{v}}_0-{\mathit{v}}_{\mathrm{grav}}$$

(48)

$${\mathit{r}}_{\mathrm{gn}}={\mathit{r}}_{\mathrm{d}}-{\mathit{r}}_0-{\mathit{r}}_{\mathrm{grav}}-\mathit{v}{t}_{\mathrm{go}}$$

(49)

where v_d and r_d are the desired terminal velocity and position in the current iteration, respectively; and ${\mathit{v}}_{\mathrm{grav}}$ and ${\mathit{r}}_{\mathrm{grav}}$ are the first and second integrals of gravitational acceleration over the time-to-go, respectively, calculated using the Taylor polynomial approximation method in spherical polar coordinates [8]. We do not elaborate upon the details of this method but instead express them as functions:

$${\mathit{v}}_{\mathrm{grav}}={\mathit{f}}_{v\mathrm{g}}\left({\mathit{v}}_0,{\mathit{v}}_{\mathrm{d}},{\mathit{r}}_0,{\mathit{r}}_{\mathrm{d}},t_{\mathrm{go}}\right)$$

(50)

$${\mathit{r}}_{\mathrm{grav}}=f_{r\mathrm{g}}\left({\mathit{v}}_0,{\mathit{v}}_{\mathrm{d}},{\mathit{r}}_0,{\mathit{r}}_{\mathrm{d}},t_{\mathrm{go}}\right)$$

(51)

By integrating the second equation in Equation system (1) twice and combining Equations (21) and (22), we obtain

$$L_{\mathrm{T}}{\mathit{\lambda }}_K+J_{\mathrm{T}}{\dot{\mathit{\lambda }}}_K={\mathit{v}}_{\mathrm{gn}}$$

(52)

$$S_{\mathrm{T}}{\mathit{\lambda }}_K+Q_{\mathrm{T}}{\dot{\mathit{\lambda }}}_K={\mathit{r}}_{\mathrm{gn}}$$

(53)

The most significant feature of explicit guidance is that only velocity constraints are used to solve the time-to-go. First, based on simplified assumptions [9]

$$\begin{array}{l}‖{\dot{\mathit{\lambda }}}_K‖\approx \omega , {\mathit{\lambda }}_K\cdot {\dot{\mathit{\lambda }}}_K\approx 0\\ \cos \left[\omega \left(t-K\right)\right]\approx 1, \sin \left[\omega \left(t-K\right)\right]\approx \omega \left(t-K\right)\end{array}$$

(54)

The simplified thrust integration coefficient can be obtained:

$$\begin{array}{l}{L}_{\mathrm{T}}\approx L_{\mathrm{P}}\\ J_{\mathrm{T}}\approx J_{\mathrm{P}}-L_{\mathrm{P}}K\end{array}$$

(55)

where

$$\begin{array}{l}{L}_{\mathrm{P}}={\displaystyle {\int }_0^{t_{\mathrm{go}}}\frac{T}{m}}\mathrm{d}t={\displaystyle {\int }_0^{t-t_{\mathrm{sum}\left(n-1\right)}}\frac{T}{m}}\mathrm{d}t+I_{s\mathrm{p}n}\ln {\displaystyle \frac{{\tau }_n}{{\tau }_n-\left(t_{\mathrm{go}}-t_{\mathrm{sum}\left(n-1\right)}\right)}}\\ J_{\mathrm{P}}={\displaystyle {\int }_0^{t_{\mathrm{go}}}\frac{T}{m}}t\mathrm{d}t\end{array}$$

(56)

The explicit expression of t_go about L_P can be obtained from Equation (56):

$$t_{\mathrm{go}}=t_{\mathrm{sum}\left(n-1\right)}+{\tau }_n-{\tau }_n{e}^{\left({\displaystyle {\int }_0^{t-t_{\mathrm{sum}\left(n-1\right)}}\frac{T}{m}}\mathrm{d}t-L_{\mathrm{P}}\right)/I_{s\mathrm{p}n}}$$

(57)

Next, letting

$$K=J_{\mathrm{P}}/L_{\mathrm{P}}$$

(58)

We obtain $J_{\mathrm{T}}\approx 0$. We can infer from Equation (52) that

$$L_{\mathrm{P}}{\mathit{\lambda }}_K\approx {\mathit{v}}_{\mathrm{gn}}$$

(59)

Considering Equation (13),

$$L_{\mathrm{P}}\approx ‖{\mathit{v}}_{\mathrm{gn}}‖$$

(60)

Then, we obtain

$$t_{\mathrm{go}}\approx t_{\mathrm{sum}\left(n-1\right)}+{\tau }_n-{\tau }_n{e}^{\left({\displaystyle {\int }_0^{t-t_{\mathrm{sum}\left(n-1\right)}}\frac{T}{m}}\mathrm{d}t-‖{\mathit{v}}_{\mathrm{gn}}‖\right)/I_{s\mathrm{p}n}}$$

(61)

Given the knowledge of v_d and r_d, ${\mathit{v}}_{\mathrm{gn}}$ is a function about t_go, so we apply numerical iteration to solve t_go from Equation (61).

Although the simplification assumptions of IGM differ from those in this paper, the final solution of t_go is obtained using Equation (61), which is the core part of its algorithm. The disadvantage of this process is that the assumptions in Formula (54) need to be valid, but these assumptions may not be met in some low-thrust and long-arc missions. Therefore, the t_go obtained using Equation (61) may also contain deviations, which is denoted as

$${\tilde{t}}_{\mathrm{go}}=t_{\mathrm{sum}\left(n-1\right)}+{\tau }_n-{\tau }_n{e}^{\left({\displaystyle {\int }_0^{t-t_{\mathrm{sum}\left(n-1\right)}}\frac{T}{m}}\mathrm{d}t-‖{\mathit{v}}_{\mathrm{gn}}‖\right)/I_{s\mathrm{p}n}}$$

(62)

In IGM, $t_{\mathrm{go}}={\tilde{t}}_{\mathrm{go}}$, which is equivalent to using the following correction formula in the main iteration of the algorithm:

$$t_{\mathrm{go}}=t_{\mathrm{go}}^{\prime }+1\cdot \left({\tilde{t}}_{\mathrm{go}}-t_{\mathrm{go}}^{\prime }\right)$$

(63)

where $t_{\mathrm{go}}^{\prime }$ represents the t_go value in the previous step of the main iteration process of the algorithm. The correction formula defaults to a step size of 1, but the algorithm may not converge due to t_go chattering in some low-thrust long-arc missions.

Unlike IGM, PEG uses the simplified ${\mathit{v}}_{\mathrm{gn}}$ in Equation (59) (usually denoted as ${\mathit{v}}_{\mathrm{go}}$) as the iteration variable for the main iteration. Therefore, t_go can be directly calculated from Equation (61) in each main iteration step. Researchers [15] modified the correction step size of ${\mathit{v}}_{\mathrm{go}}$ to 0.5 when enhancing PEG, thereby accelerating the convergence.

We used the IGM process to calculate $t_{\mathrm{go}}$. We referred to a previous study [15] for the method of enhancing PEG to eliminate chattering. We modify Equation (63) to

$$t_{\mathrm{go}}=t_{\mathrm{go}}^{\prime }+\left(1-\kappa \right)\cdot \left({\tilde{t}}_{\mathrm{go}}-t_{\mathrm{go}}^{\prime }\right)$$

(64)

where $\kappa \in \left(0, 1\right)$ is the damping factor, which was set to 0.5 in this study. The following simulation proves that this simple modification can considerably expand the adaptive boundary of the algorithm. The subiteration algorithm process for calculating t_go in DIEG is shown in

Table 1. Subiteration algorithm process for calculating tgo in DIEG.

Step	Description
1	Input $t_{\mathrm{go}}^{\prime }$ and let $t_{\mathrm{go}}^{″}=t_{\mathrm{go}}^{\prime }$.
2	Calculate ${\mathit{v}}_{\mathrm{grav}}={\mathit{f}}_{v\mathrm{g}}\left({\mathit{v}}_0,{\mathit{v}}_{\mathrm{d}},{\mathit{r}}_0,{\mathit{r}}_{\mathrm{d}},t_{\mathrm{go}}^{″}\right)$.
3	Calculate ${\mathit{v}}_{\mathrm{gn}}$ by Formula (48).
4	Calculate ${\tilde{t}}_{\mathrm{go}}$ by Formula (62).
5	Calculate $t_{\mathrm{go}}$ by Formula (64).
6	If $\left|t_{\mathrm{go}}-t_{\mathrm{go}}^{″}\right|<{\epsilon }_{\mathrm{tgo}}$, where ${\epsilon }_{\mathrm{tgo}}$ is the tolerance error, output $t_{\mathrm{go}}$; if not, make $t_{\mathrm{go}}^{″}=t_{\mathrm{go}}$ and return to Step 2.

.

4.2. Control Parameter Calculation

After obtaining the converged t_go, calculate ${\mathit{v}}_{\mathrm{grav}}$, ${\mathit{r}}_{\mathrm{grav}}$, ${\mathit{v}}_{\mathrm{gn}}$, and ${\mathit{r}}_{\mathrm{gn}}$ using Formulas (50), (51), (48), and (49), respectively. Additionally, according to the formula in ref. [6], L_P and J_P can be directly calculated, and the control parameter K can be calculated using Formula (58).

Calculate the thrust integration coefficients L_T, J_T, S_T, and Q_T using the GLFQ, and eliminate ${\dot{\mathit{\lambda }}}_K$ in Equation (52) to obtain

$${\mathit{\lambda }}_K=\mathrm{unit}\left({\mathit{v}}_{\mathrm{gn}}-{\mathit{r}}_{\mathrm{gn}}{\displaystyle \frac{J_{\mathrm{T}}}{Q_{\mathrm{T}}}}\right)$$

(65)

where ${\mathit{v}}_{\mathrm{gn}}$ and ${\mathit{r}}_{\mathrm{gn}}$ can be obtained with Formulas (48) and (49), respectively.

Before calculating ${\dot{\mathit{\lambda }}}_K$, the transversality condition given in Equation (18) must be applied to impose additional constraints on the position of the injection point. The expressions for ${\mathit{\lambda }}_{\mathrm{f}}$ and ${\dot{\mathit{\lambda }}}_{\mathrm{f}}$ can be obtained by setting $t=t_{\mathrm{go}}$ in Equation (10). In Equation (18), replace ${\mathit{r}}_{\mathrm{f}}$ and ${\mathit{v}}_{\mathrm{f}}$ with ${\mathit{r}}_{\mathrm{d}}$ and ${\mathit{v}}_{\mathrm{d}}$, respectively, to obtain

$$\mathit{A}\cdot {\mathit{\lambda }}_K=\mathit{B}\cdot {\dot{\mathit{\lambda }}}_K$$

(66)

where

$$\mathit{A}=\omega \sin \left[\omega \left(t_{\mathrm{go}}-K\right)\right]{\mathit{v}}_{\mathrm{d}}-{\displaystyle \frac{\mu \cos \left[\omega \left(t_{\mathrm{go}}-K\right)\right]}{{‖{\mathit{r}}_{\mathrm{d}}‖}^3}}{\mathit{r}}_{\mathrm{d}}$$

(67)

$$\mathit{B}=\cos \left[\omega \left(t_{\mathrm{go}}-K\right)\right]{\mathit{v}}_{\mathrm{d}}+{\displaystyle \frac{\mathsf{\mu }\sin \left[\omega \left(t_{\mathrm{go}}-K\right)\right]}{\omega {‖{\mathit{r}}_{\mathrm{d}}‖}^3}}{\mathit{r}}_{\mathrm{d}}$$

(68)

Multiply both sides of Equation (53) by $\mathit{B}$ and combine the result with Equation (66) to obtain the expression of the corrected $S_{\mathrm{T}}$:

$$S_{\mathrm{T}}={\displaystyle \frac{{\mathit{r}}_{\mathrm{gn}}\cdot \mathit{B}-Q_{\mathrm{T}}{\mathit{\lambda }}_K\cdot \mathit{A}}{{\mathit{\lambda }}_K\cdot \mathit{B}}}$$

(69)

Finally, the expression for ${\dot{\mathit{\lambda }}}_K$ is obtained:

$${\dot{\mathit{\lambda }}}_K={\displaystyle \frac{{\mathit{r}}_{\mathrm{gn}}-S_{\mathrm{T}}{\mathit{\lambda }}_K}{Q_{\mathrm{T}}}}$$

(70)

After both ${\mathit{\lambda }}_K$ and ${\dot{\mathit{\lambda }}}_K$ are calculated, calculate the values of D₁ and D₂ according to Formula (24) and use the GLFQ again to calculate the thrust integration coefficients L_T, J_T, S_T, and Q_T after correcting D₁ and D₂.

4.3. Prediction and Correction of Injection Points

${\mathit{v}}_{\mathrm{thrust}}$ and ${\mathit{r}}_{\mathrm{thrust}}$ can be calculated according to Formulas (21) and (22), and the predicted terminal velocity ${\mathit{v}}_{\mathrm{p}}$ and position ${\mathit{r}}_{\mathrm{p}}$ are

$$\begin{array}{l}{\mathit{v}}_{\mathrm{p}}={\mathit{v}}_0+{\mathit{v}}_{\mathrm{grav}}+{\mathit{v}}_{\mathrm{thrust}}\\ {\mathit{r}}_{\mathrm{p}}={\mathit{r}}_0+{\mathit{v}}_0{t}_{\mathrm{go}}+{\mathit{r}}_{\mathrm{grav}}+{\mathit{r}}_{\mathrm{thrust}}\end{array}$$

(71)

Project ${\mathit{r}}_{\mathrm{p}}$ onto the target orbital plane, and the resulting vector is denoted as ${\overline{\mathit{r}}}_{\mathrm{p}}$. Let the intersection point between ${\overline{\mathit{r}}}_{\mathrm{p}}$ and the target orbit be M; then, M is the target injection point for the next iteration step, that is, ${\mathit{r}}_{\mathrm{d}}=\overrightarrow{OM}$, as shown in Figure 2.

Let the angle between ${\mathit{r}}_{\mathrm{d}}$ and the OY axis be $\theta$, and its expression is

$$\theta =\arctan \left({\overline{r}}_{\mathrm{p}X}/{\overline{r}}_{\mathrm{p}Y}\right)$$

(72)

The true anomaly of the target injection point is

$$f_{\mathrm{d}}=f^{\ast }+\theta$$

(73)

The modulus of ${\mathit{r}}_{\mathrm{d}}$ is

$$‖{\mathit{r}}_{\mathrm{d}}‖={\displaystyle \frac{a^{\ast }\left(1-e^{\ast 2}\right)}{1+e^{\ast }\cos f_{\mathrm{d}}}}$$

(74)

Therefore,

$${\mathit{r}}_{\mathrm{d}}=\left[\begin{array}{l}{r}_{\mathrm{d}x}\\ r_{\mathrm{d}y}\\ r_{\mathrm{d}z}\end{array}\right]=\left[\begin{array}{c}‖{\mathit{r}}_{\mathrm{d}}‖\sin \theta \\ ‖{\mathit{r}}_{\mathrm{d}}‖\cos \theta \\ 0\end{array}\right]$$

(75)

Let the vertical velocity at point M relative to the ground be ${\overline{v}}_{\mathrm{d}y}$ and the horizontal velocity be ${\overline{v}}_{\mathrm{d}x}$. Their expressions are

$$\begin{array}{l}{\overline{v}}_{\mathrm{d}x}=\sqrt{\mu a^{\ast }\left(1-e^{\ast 2}\right)}/‖{\mathit{r}}_{\mathrm{d}}‖\\ {\overline{v}}_{\mathrm{d}y}=e^{\ast }\sin f_{\mathrm{d}}\sqrt{\mu /a^{\ast }\left(1-e^{\ast 2}\right)}\end{array}$$

(76)

We obtain ${\mathit{v}}_{\mathrm{d}}$ for the next main iteration step:

$${\mathit{v}}_{\mathrm{d}}=\left[\begin{array}{c}{\overline{v}}_{\mathrm{d}x}\cos \theta +{\overline{v}}_{\mathrm{d}y}\sin \theta \\ -{\overline{v}}_{\mathrm{d}x}\sin \theta +{\overline{v}}_{\mathrm{d}y}\cos \theta \\ 0\end{array}\right]$$

(77)

During initialization, let $\theta =0$ to obtain ${\mathit{r}}_{\mathrm{d}}$ and ${\mathit{v}}_{\mathrm{d}}$ for the first iteration. Let the $r_{\mathrm{d}x}$ obtained from the previous iteration be denoted as $r_{\mathrm{d}x}$, and $r_{\mathrm{d}x}^{\prime }=0$ during initialization. When $\left|r_{\mathrm{d}x}-r_{\mathrm{d}x}^{\prime }\right|\le {\epsilon }_{r\mathrm{d}x}$, DIEG has converged, where ${\epsilon }_{r\mathrm{d}x}$ is the tolerance error. The current guidance command ${\mathit{u}}_0$ is given in Equation (78):

$${\mathit{u}}_0=\mathrm{unit}\left({\displaystyle \frac{{\mathit{\lambda }}_K\cos \left(\omega K\right)-\frac{{\dot{\mathit{\lambda }}}_K}{\omega }\sin \left(\omega K\right)}{‖{\mathit{\lambda }}_K\cos \left(\omega K\right)-\frac{{\dot{\mathit{\lambda }}}_K}{\omega }\sin \left(\omega K\right)‖}}\right)$$

(78)

The pitch angle command ${\varphi }_G$ and yaw angle command ${\psi }_G$ in the launch inertial system are given in Equation (79):

$$\begin{array}{l}{\mathit{u}}_{0\mathrm{G}}=\left[\begin{array}{l}{u}_{0\mathrm{G}x}\\ u_{0\mathrm{G}y}\\ u_{0\mathrm{G}z}\end{array}\right]={\mathit{M}}_{\mathrm{GO}}\cdot {\mathit{u}}_0\\ {\varphi }_G=\arctan \left({\displaystyle \frac{u_{0\mathrm{G}y}}{u_{0\mathrm{G}x}}}\right)\\ {\psi }_G=\arcsin \left(-u_{0\mathrm{G}z}\right)\end{array}$$

(79)

where ${\mathit{M}}_{\mathrm{GO}}$ is the transformation matrix from the reference injection point coordinate system to the launch inertial coordinate system.

4.4. Iterative Process of DIEG

The complete iterative process of DIEG is shown in

Table 2. Iterative process of DIEG.

Step	Description
1	Initialize $t_{\mathrm{go}}^{\prime }$, ${\mathit{r}}_{\mathrm{d}}$, ${\mathit{v}}_{\mathrm{d}}$, $r_{\mathrm{d}x}^{\prime }$, $D_1$, and $D_2$.
2	At the beginning of a new flight phase, input parameters for each phase $I_{\mathrm{sp}i}$, ${\tau }_i$, and $t_{\mathrm{b}i}$; input the remaining number of phases n.
3	At each guidance period, input the current flight states ${\mathit{r}}_0$, ${\mathit{v}}_0$, and $a_{\mathrm{c}0}$.
4	Calculate $t_{\mathrm{go}}$ by the method described in Section 4.1.
5	Calculate Lp and Jp using the previously described method [6], and calculate $K$ using Formula (58).
6	Calculate LT, JT, ST, and QT by using the GLQF method described in Section 3.
7	Calculate ${\mathit{\lambda }}_K,{\dot{\mathit{\lambda }}}_K$ using the method described in Section 4.2.
8	Update $D_1$ and $D_2$ using Formula (24), and calculate LT, JT, ST, and QT again.
9	Predict and correct the position of the injection point by the method described in Section 4.3.
10	If $\left|r_{\mathrm{d}x}-r_{\mathrm{d}x}^{\prime }\right|\le {\epsilon }_{r\mathrm{d}x}$, calculate ${\varphi }_G$ and ${\psi }_G$ using Formula (79); if not, return to Step 4.

. DIEG is not used when $t_{\mathrm{go}}\le 5 \mathrm{s}$ to prevent algorithm divergence due to insufficient position error in the last few seconds of flight, and the attitude angle remains unchanged until the semimajor axis satisfies the constraint, at which point the engine shuts down.

5. Simulation and Analysis

This section describes the simulation of the three phases of a rocket before entering the target orbit. The rocket parameters for each flight phase are shown in

Table 3. Rocket parameters for each phase.

	Phase 1	Phase 2	Phase 3
Mass of fuel (kg)	168,485.0	58,252.0	58,252.0
Mass of rocket structure (kg)	116,844.0	38,743.0	38,743.0
Magnitude of thrust (N)	2,869,956	0	271,207
Specific impulse (m/s)	3389.7	0	4410.6
Flight time (s)	199.00	11.00	263.45

, with Phase 2 being a coasting phase. The orbital elements of the initial point and nominal injection point of the simulation are shown in

Table 4. Orbital elements of initial and nominal injection points.

Orbital Element	Initial Point in Multiphase Mission *	Initial Point in Single-Phase Mission **	Nominal Injection Point
a (m)	3,774,423.02	5,500,908.72	6,560,590.37
e	0.734989	0.195333	0.002669
i (deg)	24.59149	26.61440	26.95575
Ω (deg)	342.53390	338.01405	337.61593
w (deg)	313.70517	327.20465	290.19344
f (deg)	176.21225	176.66212	231.51759

* The initial point in Phase 1; ** the initial point in Phase 3.

. The meanings of the five algorithms mentioned in the following text are as follows:

• hp-RPM: hp adaptive Radau pseudo-spectral method, implemented through GPOPS [41], which is used as the optimal algorithm in this study.
• Algorithm 0: The DIEG with ${\epsilon }_{\mathrm{tgo}}$ and ${\epsilon }_{r\mathrm{d}x}$ set to 0.01 s and 100 m, respectively.
• Algorithm 1: Set damping factor $\kappa =0$ in DIEG.
• Algorithm 2: Replace GLQF with linearized thrust integral [7] in DIEG.
• Algorithm 3: Replace the transversality condition (18) with the orthogonal assumption (19) in DIGM.
• Algorithm 4: PEG in ref. [7], but the gravity integration method in ref. [8] was adopted, similar to the other algorithms mentioned above, and the update step size of ${\mathit{v}}_{\mathrm{go}}$ was set to 0.5, as in reference [15].

5.1. Simulation Results and Analysis of Single-Phase Mission

The above five methods were used for numerical simulation starting from 100% and gradually decreasing the engine thrust at intervals of one, and five representative scenarios were obtained:

• 2%: the maximum thrust drop ratio where Algorithm 4 converges;
• 9%: the maximum thrust drop ratio where Algorithm 1 converges;
• 14%: the maximum thrust drop ratio where Algorithm 2 converges;
• 34%: the maximum thrust drop ratio where Algorithm 3 converges;
• 35%: the maximum thrust drop ratio where Algorithm 0 converges.

The comparison of the command angle among the above five scenarios is shown in Figure 3, Figure 4, Figure 5, Figure 6 and Figure 7; the terminal accuracy, convergence speed, and flight time are compared in

Table 5. Performance of different algorithms in various thrust reduction ratio cases in the single-phase mission.

Thrust Drop Ratio	Algorithm Number	Terminal Deviation				Maximum Number of Iteration Steps	Flight Time Longer Than Optimal Solution (s)
		e	i (deg)	Ω (deg)	w (deg)
2%	0	−3.36 × 10−6	1.67 × 10−5	−1.24 × 10−5	0.037	7	0.002
	1	−3.35 × 10−6	1.67 × 10−5	−1.24 × 10−5	0.036	11	0.002
	2	−8.10 × 10−6	1.05 × 10−4	−7.77 × 10−5	0.128	7	0.015
	3	−2.61 × 10−6	1.45 × 10−5	−1.07 × 10−5	0.038	7	0.001
	4	−2.65 × 10−6	1.38 × 10−5	−1.02 × 10−5	0.038	6	0.001
9%	0	−3.87 × 10−6	4.27 × 10−5	−2.92 × 10−5	0.045	9	0.002
	1	−3.86 × 10−6	4.26 × 10−5	−2.91 × 10−5	0.044	560	0.002
	2	−9.32 × 10−6	1.23 × 10−4	−8.37 × 10−5	0.149	8	0.029
	3	2.89 × 10−6	3.69 × 10−5	−2.52 × 10−5	0.045	9	0.001
14%	0	−4.03 × 10−6	5.72 × 10−5	−3.62 × 10−5	0.055	11	0.003
	2	−1.04 × 10−5	1.29 × 10−4	−8.18 × 10−5	0.158	7	0.052
	3	−3.48 × 10−6	4.92 × 10−5	−3.12 × 10−5	0.039	11	0.004
34%	0	−5.70 × 10−6	8.91 × 10−5	−3.13 × 10−5	0.067	17	0.002
	3	−4.00 × 10−6	6.37 × 10−5	−2.24 × 10−5	0.036	114	0.155
35%	0	−6.09 × 10−6	8.97 × 10−5	−2.98 × 10−5	0.059	59	0.002

. Algorithm 0 has the largest convergence range, and its results are almost the same as those of hp-RPM, proving its strong convergence and optimality. Algorithm 1 has similar accuracy and optimality to Algorithm 0, but the convergence range and convergence speed are substantially reduced, indicating that introducing the damping factor to the iterative calculation of t_go considerably improves the convergence of the algorithm. The accuracy and optimality of Algorithm 2 are the lowest, and the convergence range is smaller than that of Algorithm 0, which verifies the necessity of introducing an accurate thrust integration method. The performance of Algorithm 3 is slightly better than that of Algorithm 0 when the degree of failure is small, which indicates that the orthogonal assumption used in the previous PEG algorithm is rational. However, the optimality gradually weakens compared with that of Algorithm 0 as the thrust further decreases, and the range of faults to which it can adapt is also smaller than that of Algorithm 0. The adaptability range of Algorithm 4 is the smallest, indicating the limitations of the classical PEG algorithm in dealing with fault conditions, also reflecting the need to improve explicit guidance.

Figure 8 compares the t_go curves of the above algorithm during the first guidance period when the thrust drop ratio is 34%. Algorithm 0 quickly converged without notable chattering. The result of Algorithm 1 oscillated with equal amplitude without any convergence trend, which was caused by a lack of damping. After several tens of large oscillations, the result of Algorithm 2 oscillated slightly near the value with a substantial error from the optimal value, indicating that low-accuracy linearized thrust integration results in a large error in the solution, ultimately leading to nonconvergence. The results of Algorithm 3 were almost consistent with those of Algorithm 1 before Step 8 but then widely oscillated and gradually converged, indicating that the accuracy of the transversality condition weakly impacted the performance of the algorithm, but this impact increased at the convergent boundary. The result of Algorithm 4 slowly converged to an error value that was close to that of Algorithm 2, which was caused by inaccurate linearization integration. Although Algorithm 4 converged in the first step, divergence still occurred later due to excessive deviation from the optimal solution.

5.2. Simulation Results and Analysis of Multiphase Mission

The above five methods were used for numerical simulation in the multiphase mission, starting from 100% and gradually decreasing the engine thrust at intervals of 1%. Five scenarios were obtained:

• 1%: the maximum thrust drop ratio where Algorithm 4 converges;
• 11%: the maximum thrust drop ratio where Algorithm 1 converges;
• 49%: the maximum thrust drop ratio where Algorithm 3 converges;
• 51%: the maximum thrust drop ratio where Algorithm 2 converges;
• 55%: the maximum thrust drop ratio where Algorithm 0 converges.

The command angle in the above five scenarios is compared in Figure 9, Figure 10, Figure 11, Figure 12 and Figure 13, and the terminal accuracy, convergence speed, and flight time are compared in

Table 6. Performance of different algorithms under various thrust reduction ratios in the multiphase mission.

Thrust Drop Ratio	Algorithm Number	Terminal Deviation				Maximum Number of Iteration Steps	Flight Time Longer Than Optimal Solution (s)
		e	i (deg)	Ω (deg)	w (deg)
1%	0	−8.67 × 10−7	−3.11 × 10−5	2.33 × 10−5	−0.001	9	0.024
	1	−9.83 × 10−6	1.85 × 10−4	1.44 × 10−4	0.170	15	0.030
	2	−6.18 × 10−6	6.94 × 10−5	−5.24 × 10−5	0.107	9	0.109
	3	−2.87 × 10−7	4.57 × 10−5	−3.42 × 10−5	−0.019	9	0.002
	4	−2.30 × 10−6	4.74 × 10−6	−3.65 × 10−6	0.030	146	0.014
11%	0	−9.76 × 10−7	5.32 × 10−6	3.81 × 10−6	−0.003	10	0.042
	1	−1.21 × 10−5	1.65 × 10−4	1.19 × 10−4	0.189	118	0.049
	2	−8.19 × 10−6	8.69 × 10−5	−6.05 × 10−5	0.131	10	0.215
	3	8.48 × 10−7	−1.40 × 10−5	9.66 × 10−6	−0.018	10	0.002
49%	0	−6.20 × 10−6	1.75 × 10−4	−3.16 × 10−5	−0.046	21	0.374
	2	3.32 × 10−5	2.72 × 10−4	−4.06 × 10−5	−0.308	31	10.525
	3	3.62 × 10−6	1.51 × 10−4	−2.78 × 10−5	−0.039	21	0.572
51%	0	−7.14 × 10−6	1.87 × 10−4	−2.16 × 10−5	0.052	21	0.384
	2	−4.39 × 10−5	3.99 × 10−4	−2.87 × 10−5	−0.373	22	14.849
55%	0	−9.81 × 10−6	2.18 × 10−4	1.17 × 10−5	0.065	73	0.273

. Algorithms 0 and 4 have the strongest and weakest adaptability, respectively, further verifying that the algorithm proposed in this paper can increase the adaptability of rockets to thrust faults. Although Algorithm 2 adapts to a range of faults, similar to Algorithm 0, its deviation from the optimal solution is the largest. From the perspectives of terminal accuracy and optimality, Algorithm 3 performs the best under small faults, but command angle chatters as in Algorithm 1 when approaching the injection point, which poses a risk in their application. Although the optimality of DIEG decreases with increasing fault severity in the multiphase mission, the results of the simulation experiments showed that when the thrust drop ratio does not exceed 55%, the difference in flight time between DIEG and the optimal solution is less than 0.384 s, indicating strong optimality.

We conducted 1000 Monte Carlo shooting simulations to verify the robustness of DIEG to the various parameter deviations shown in

Table 7. Setting of error terms in Monte Carlo simulation.

Deviation Term	3σ Value in Phase 1	3σ Value in Phase 3
Initial position (km)	(20, 10, 5)	-
Initial velocity (m/s)	(100, 50, 20)	-
Specific impulse (m/s)	15	15
Deviation ratio of fuel mass flow rate	3%	3%
Deviation ratio of fuel mass	0.5%	0.5%
Deviation ratio of structure mass	1%	1%

. The simulation results are shown in Figure 14, Figure 15 and Figure 16 and

Table 8. Statistics of orbit injection accuracy.

Deviation Term of Orbital Elements	Mean	Standard Deviation
a (m)	2.821	1.578
e	1.406 × 10−6	3.310 × 10−6
i (deg)	7.223 × 10−5	7.138 × 10−5
Ω (deg)	5.449 × 10−5	5.468 × 10−5
w (deg)	−3.558 × 10−2	−5.712 × 10−2
f (deg)	7.378 × 10−4	0.6983

. Figure 14 shows that the command angle suddenly changes at the beginning of Phase 3 because of the deviation in the fuel mass and its flow rate in Phase 1, which causes the engine to shut down before the preset t_b1. However, the command angle curve remains smooth during each phase. The statistical results in Figure 14 and Table 8 indicate that using DIEG enables the rocket to accurately enter orbit with five terminal constraints under various deviations. Figure 16 shows that the maximum iteration steps of the algorithm are mostly distributed within the range of [8,11], with a maximum value of 12, indicating that the algorithm still quickly converges under parameter deviation. Only 1–2 steps are required to converge in the remaining guidance periods after the algorithm converges in the first guidance period of each burn phase. The above results collectively demonstrate that DIEG is robust.

6. Conclusions

We developed a DIEG algorithm to increase the fault tolerance of multistage rockets to thrust drop faults. In the algorithm, we replaced the simplified assumptions with the accurate transversality condition, calculated thrust integration using GLQF instead of linearization methods, and introduced damping factors when solving time-to-go. The following conclusions were drawn by comparing DIEG with PEG and an optimization algorithm, as well as conducting Monte Carlo simulations under deviations:

(1) DIEG increases the maximum thrust drop ratio to which the guidance algorithm can adapt from 2% for PEG to 35% and from 1% to 55% in a single- and multiphase mission, respectively, which shows its stronger convergence.

(2) With DIEG, the rocket flies no more than 3 ms longer than the optimal solution on a single-phase mission; DIEG allows the rocket to fly no more than 0.384 s longer than the optimal solution when the thrust drops by no more than 55% in a multiphase mission, demonstrating its strong optimality.

(3) The DIEG algorithm can converge within 12 steps under various deviations, and the terminal accuracy meets the requirements, indicating its robustness and practical engineering value.

(4) Among the three improvement measures applied in this study, the introduction of the damping factor most strongly impacted algorithm performance, which enhances the convergence of the algorithm by eliminating the chattering in time-to-go. The second strongest impact is the improvement in thrust integration, which enhances the algorithm’s optimality by increasing the accuracy of prediction and correction, thereby strengthening the algorithm convergence. The weakest impact was achieved with the improvement in the transversality condition, which increasingly impacts the algorithm as the degree of failure increases.