Home Journals Articles Updates About
CN
image
收藏切换
Online reentry guidance algorithm based on trajectory analytical solutions
收藏切换
PDF
Weibo SUN1, 2, Ping MA1, 2, Xiaonan LI1, 2, Songyan WANG1, 2, Tao CHAO1, 2, *
Journal of Systems Engineering and Electronics | 2026, 37(3) : 1002 - 1018
Less
收藏切换
Journal of Systems Engineering and Electronics | 2026, 37(3): 1002-1018
• CONTROL THEORY AND APPLICATION •
Online reentry guidance algorithm based on trajectory analytical solutions
Full
Weibo SUN1, 2, Ping MA1, 2, Xiaonan LI1, 2, Songyan WANG1, 2, Tao CHAO1, 2, *
Affiliations
  • 1Control and Simulation Center, Harbin Institute of Technology, Harbin 150001, China
  • 2National Key Laboratory of Modeling and Simulation for Complex Systems, Harbin 150001, China
Published: 2026-06-18 doi: 10.23919/JSEE.2026.000089
Outline
收藏切换

To enhance the real-time performance and accuracy of guidance command generation, we propose an online reentry guidance algorithm based on analytical solutions of the hypersonic glide trajectory (HGT). Initially, an altitude-velocity profile is designed in the longitudinal plane to satisfy both path and terminal constraints. Based on this profile, we derive analytical solutions for the flight path angle (FPA) and bank angle. Subsequently, by employing the Newton-Raphson method to linearize the reentry motion equations, analytical solutions for the latitude and heading angle are obtained. Furthermore, we introduce an improved particle swarm optimization (IPSO) algorithm to optimize the profile parameters. This approach significantly enhances the algorithm’s global convergence by narrowing the parameter optimization range and adaptively adjusting the inertia weight and cognitive factors. Finally, we present an online guidance algorithm that combines the HGT analytical solutions with the IPSO algorithm. This algorithm effectively achieves longitudinal and lateral guidance by continuously updating the altitude-velocity profile and bank angle symbol in real time. Simulation results demonstrate that the proposed algorithm is fast, efficient, accurate, and holds significant potential for broader application.

altitude-velocity profile  /  analytical solution  /  improved particle swarm optimization (IPSO)  /  online reentry guidance
Weibo SUN, Ping MA, Xiaonan LI, Songyan WANG, Tao CHAO. Online reentry guidance algorithm based on trajectory analytical solutions[J]. Journal of Systems Engineering and Electronics, 2026 , 37 (3) : 1002 -1018 . DOI: 10.23919/JSEE.2026.000089
Full Text
Less
1.　Introduction
Less
Hypersonic glide vehicles (HGVs) are capable of performing high-speed flight over a wide range due to their unique aerodynamic design. These vehicles have a wide range of applications, including military reconnaissance, rapid strike, space flight, and even potential future civilian high-speed transportation [1,2]. They usually fly at the edge of the atmosphere or in higher airspace and can flexibly adjust their flight paths in response to changes in atmospheric density. This ability enables them to execute rapid maneuvers and directional changes. Such flight flexibility places extremely high technical demands on the design of the vehicles’ reentry guidance methods.
Reentry guidance algorithms are predominantly classified into two categories: standard trajectory guidance (STG) and predictor-corrector guidance (PCG) [3,4]. STG, the original algorithm developed, involves offline trajectory design and online trajectory tracking. Standard trajectories are precomputed and stored in the onboard computer using planning or optimization techniques. A stable tracking controller is then designed to handle disturbances and uncertainties. Over recent decades, STG and its improvements have been extensively employed in various reentry vehicles [5–9]. He et al. [5] effectively reduced computational demands by employing the Newton iteration method for calculating longitudinal and lateral profiles. Additionally, they developed a proportional integral derivative (PID) controller that exhibits greater robustness compared to traditional tracking methods. Morio et al. [6] introduced an innovative approach for generating standard trajectories in the longitudinal plane, utilizing the flatness approach to mitigate disturbances through the PID controller. Yan et al. [7] designed a state-feedback law that incorporates drag and its rate of change as feedback to track the standard drag profile. Mease et al. [8,9] proposed a trajectory planner that generates feasible trajectories and corresponding bank angle profiles. Despite its widespread use, STG, which relies on precomputed standard trajectories, is unable to adapt to real-time environmental changes, rendering it unsuitable for online guidance systems that require real-time responsiveness to external inputs.
PCG adjusts the guidance commands based on the deviation between the terminal position of the actual trajectory and that of the predicted trajectory. Wang et al. [10] introduced a novel method utilizing a piecewise linear bank profile design, which not only strictly constrains the terminal altitude but also ensures adherence to range constraints, effectively addressing the mismatch between terminal range and altitude. Additionally, Wang et al. [11,12] proposed a method that incorporates no-fly zone and waypoint constraints. By employing numerical integration and fuzzy logic, this method effectively predicts and corrects the trajectory, demonstrating enhanced robustness and operational efficiency. Zhang et al. [13] developed an improved PCG method that significantly enhances the guidance accuracy and robustness of high lift-to-drag ratio reusable vehicles through parametric design and feedback correction. Despite its adaptability to dynamic environmental disturbances and uncertainties, the substantial computational demands of PCG may limit its applicability in vehicles with constrained computational resources.
Recently, several novel online guidance algorithms have been developed. Xu et al. [14] proposed a reentry guidance framework that significantly improved real-time performance, accuracy, and robustness of guidance command generation. This was achieved through parameterized processing and swarm intelligence optimization algorithms. Pan et al. [15] introduced a three-dimensional guidance method by constructing a precise entry guidance model and a multi-constraint flight corridor, which enhanced the vehicle’s capability for long-distance strikes and extensive lateral maneuvers. Yu et al. [16–18] developed a computationally efficient analytical reentry guidance method that rapidly plans the trajectory using improved analytical solutions and addresses coupling issues within motion equations. Building on Yu’s work, Yang et al. [19] applied a Newton-like method to solve highly nonlinear models of the reentry motion equations and established a Chebyshev series analytical formula for the rapid and efficient prediction of three-dimensional hypersonic gliding trajectories. Zeng et al. [20] introduced a new three-dimensional guidance algorithm that reduces computational load and simplifies attitude control by optimizing bank angle reversals and streamlining the motion equations. Zhou et al. [21,22] developed a new glide guidance law that quickly generates guidance commands by designing separate longitudinal and lateral paths, ensuring compliance with terminal altitude, flight path angle (FPA), and position requirements. By incorporating virtual targets and applying analytical solutions, they simplified the implementation of the guidance law, enhancing its robustness. Additionally, Zhou et al. [23] proposed a novel lateral guidance strategy that integrates analytical solutions with reinforcement learning algorithms, effectively addressing complex trajectory constraints and uncertain disturbances while optimizing performance and accuracy. Despite these advancements, the new methods still face challenges in meeting the rigorous demands of practical engineering applications.
In order to address these issues, we propose an online guidance algorithm based on analytical solutions of the hypersonic glide trajectory (HGT). Initially, a fifth-order polynomial altitude-velocity profile is designed. Based on this profile, we derive analytical solutions for the FPA and bank angle. Subsequently, we simplify the reentry motion equations and apply the Newton-Raphson method to linearize these simplified equations, yielding analytical solutions for the latitude and heading angle. Next, we determine the profile parameters based on initial and terminal conditions on altitude, velocity, and FPA, as well as two specially designed altitude parameters. Concurrently, we make improvements to the particle swarm optimization (PSO) algorithm. These improvements include narrowing the parameter optimization range and adaptively adjusting inertial weights and cognitive factors, which significantly accelerate the global convergence of the algorithm. Improved PSO (IPSO) algorithm, which takes the terminal range as the objective function and considers path constraints, is used to optimize two altitude parameters to determine the profile. Finally, we propose an online guidance algorithm based on the HGT analytical solutions and the IPSO algorithm and we achieve online guidance of the vehicle by optimizing the altitude-velocity profile for each guidance cycle and combining it with the bank reversal strategy. The main contributions of this paper are summarized as follows:
(i) The HGT analytical solutions are derived, which simultaneously solves the longitudinal and lateral motion states analytically, significantly reducing the reentry motion equations to first order.
i) Analytical formulas for the FPA and bank angle, with the altitude-velocity profile as key parameters, are derived to accommodate various trajectory forms and fully consider the influence of the Earth’s rotation to ensure trajectory accuracy.
ii) Based on the Newton-Raphson method, the lateral motion equations are linearized and transformed into the linear recursive equations, from which the recursive analytical formulas for latitude and heading angle are derived.
(ii) To achieve online guidance, a simplified profile updating strategy is proposed to replace the traditional trajectory tracker design, along with an online guidance algorithm based on the IPSO algorithm.
i) A fifth-order polynomial altitude-velocity profile is designed, with coefficients determined by two altitude parameters, which naturally satisfy the constraints of terminal altitude, velocity, and FPA. The vehicle is controlled by optimizing these altitude parameters to meet the terminal range requirements.
ⅱ) Path constraints are introduced to shorten the particle position boundary. To address the issue of frequent particle movement when approaching the optimal solution, a strategy to reduce particle speed is adopted, allowing the speed boundary to shrink linearly with iteration time, thereby improving search speed.
ⅲ) An adaptive update strategy for inertia weight and cognitive factors is designed to further improve search efficiency.
The remainder of this paper is structured as follows: Section 2 describes the reentry guidance problems. Section 3 introduces the derivation process of the HGT analytical solutions. Section 4 elaborates on the framework of the online reentry guidance algorithm, including the profile updating strategy and the IPSO. Section 5 verifies the performance of the reentry guidance algorithm through various simulation examples. Section 6 concludes this paper and discusses future research directions.
2.　Problem formulation
Less
2.1　Reentry motion equations
Considering the HGV as a point particle and the Earth as a uniform sphere, the reentry motion equations of the HGV can be described as
$ \frac{\mathrm{d}H}{\mathrm{d}t}=V\sin\gamma, $
$ \frac{\mathrm{d}\lambda}{\mathrm{d}t}=\frac{V\cos\gamma\sin\psi}{R\cos\phi}, $
$ \frac{\mathrm{d}\phi}{\mathrm{d}t}=\frac{V\cos\gamma\cos\psi}{R}, $
$ \begin{split}& \frac{\mathrm{d}V}{\mathrm{d}t}=-\frac{D}{m}-g\sin\gamma+R\omega_{\text{e}}^2\cos\phi\cdot \\& \; \; \left(\sin\gamma\cos\phi-\cos\psi\cos\gamma\sin\phi\right),\end{split} $
$\begin{split}& \frac{\mathrm{d}\gamma}{\mathrm{d}t}=\frac{L\cos\sigma}{mV}+\left(\frac{V}{R}-\frac{g}{V}\right)\cos\gamma+2\omega_{\text{e}}\cos\phi\sin\psi+ \\& \qquad\frac{R\omega_{\text{e}}^2}{V}\left(\cos\gamma\cos^2\phi+\cos\psi\sin\gamma\sin\phi\cos\phi\right), \end{split}$
$ \begin{split} & \frac{\mathrm{d}\psi}{\mathrm{d}t}=\frac{L\sin\sigma}{mV\cos\gamma}+2\omega_{\text{e}}\sin\phi+\frac{V}{R}\sin\psi\tan\phi\cos\gamma+\; \\ & \quad\frac{R\omega_{\text{e}}^2\sin\psi\sin\phi\cos\phi}{V\cos\gamma}-\frac{2\omega_{\text{e}}\cos\psi\cos\phi\sin\gamma}{\cos\gamma},\end{split} $
where $R$ is the radial distance, $R = {R_{\text{e}}} + H$, ${R_{\text{e}}}$ is the Earth radius, $H$ is the altitude, $\lambda $ is the longitude, $\phi $ is the latitude, $V$ is the velocity, $\gamma $ is the FPA, $\psi $ is the heading angle, $g$ is the acceleration of gravity, $\sigma $ is the bank angle, $m$ is the mass, ${\omega _{\text{e}}}$ is the angular rotation velocity of the Earth, $L$ is the aerodynamic lift force, and $D$ is the aerodynamic drag force.
$L$ and $D$ are as follows:
$ \left\{\begin{gathered}L=0.5\rho V^2SC_L\left(\alpha,\mathrm{Ma}\right) \\ D=0.5\rho V^2SC_D\left(\alpha,\mathrm{Ma}\right) \\ \end{gathered}\right. $
where $S$ is the reference area, $\rho $ is the atmospheric density, ${C_L}$ is the aerodynamic lift coefficient, ${C_D}$ is the aerodynamic drag coefficient, both coefficients are functions of the angle of attack (AOA) $\alpha $ and Mach ${\mathrm{Ma}}$ as independent variables, ${\mathrm{Ma}} = {V \mathord{\left/ {\vphantom {V {{V_{\text{s}}}}}} \right. } {{V_{\text{s}}}}}$, and ${V_{\text{s}}}$ is the sound velocity.
The atmosphere model is
$ \rho = {\rho _0} \cdot \exp \left( { - \beta H} \right) $
where $ {\rho _0} = 1.225\;{\text{kg/}}{{\text{m}}^3} $ and $ \beta = 1.406 \times {10^{ - 4}}\;{{\text{m}}^{ - 1}} $.
2.2　Reentry constraints
When the HGV reenters the atmosphere, its flight trajectory is subject to the following path constraints:
$ \left\{ \begin{gathered} \dot Q = {K_{\dot Q}}{\rho ^{0.5}}{V^{3.15}} \leqslant {{\dot Q}_{\max }} \\ n = {{\sqrt {{L^2} + {D^2}} } \mathord{\left/ {\vphantom {{\sqrt {{L^2} + {D^2}} } {mg}}} \right. } {(mg)}} \leqslant {n_{\max }} \\ q = 0.5\rho {V^2} \leqslant {q_{\max }} \\ \end{gathered} \right. $
where $\dot Q$ represents the heat flow, $n$ represents the overload, $q$ represents the dynamic pressure, ${\dot Q_{\max }}$, ${n_{\max }}$, and ${q_{\max }}$ are the maximum allowable values of the path constraints. $ {K_{\dot Q}} $ is the coefficient of the heat flow model, determined by the aerodynamic structure of the HGV.
The terminal constraints of the HGV are
$ \left\{\begin{aligned}&H\left( {{t_{\text{f}}}} \right) = {H_{\text{f}}}\\&V\left( {{t_{\text{f}}}} \right) = {V_{\text{f}}}\\&\gamma \left( {{t_{\text{f}}}} \right) = {\gamma _{\text{f}}}\\&{R_{{\text{stogo}}}}\left( {{t_{\text{f}}}} \right) \leqslant {R_{{\text{stogof}}}}\end{aligned}\right. $
where $ {H_{\text{f}}} $, $ {V_{\text{f}}} $, $ {\gamma _{\text{f}}} $ and $ {R_{{\text{stogof}}}} $ represent the expected terminal values at the terminal time $ {t_{\text{f}}} $, $ {R_{{\text{stogo}}}} $ is the range-to-go.
$ {R_{{\text{stogo}}}} = {R_{\text{e}}}\arccos \left[ {\sin {\phi _{\text{T}}}\sin \phi + \cos {\phi _{\text{T}}}\cos \phi \cos \left( {{\lambda _{\text{T}}} - \lambda } \right)} \right] $
where $ {\lambda _{\text{T}}} $ and $ {\phi _{\text{T}}} $ are the longitude and latitude of target.
In addition, the control constraints should also be satisfied.
$ \left\{\begin{aligned}&\alpha \leqslant {\alpha _{\max }}\\&- {\sigma _{\max }} \leqslant \sigma \leqslant {\sigma _{\max }}\end{aligned}\right. $
where $ {\alpha _{\max }} $ and $ {\sigma _{\max }} $ are the maximum allowable values of the control variables.
3.　HGT analytical solutions
Less
In this paper, we transform the reentry motion equations into a form with $V$ as the independent variable.
$ \frac{{{\mathrm{d}}H}}{{{\mathrm{d}}V}} = \frac{{{\mathrm{d}}t}}{{{\mathrm{d}}V}}V\sin \gamma $
$\begin{split}& \frac{{{\mathrm{d}}\gamma }}{{{\mathrm{d}}V}} = \frac{{{\mathrm{d}}t}}{{{\mathrm{d}}V}}\Bigg[ {\frac{{L\cos \sigma }}{{mV}} + \left( {\frac{V}{R} - \frac{g}{V}} \right)\cos \gamma + 2{\omega _{\text{e}}}\cos \phi \sin \psi + } \\&\qquad {\frac{{R\omega _{\text{e}}^2\left( {\cos \gamma {{\cos }^2}\phi + \cos \psi \sin \gamma \sin \phi \cos \phi } \right)}}{V}} \Bigg]\\[-1pt] \end{split} $
$ \frac{{{\mathrm{d}}\lambda }}{{{\mathrm{d}}V}} = \frac{{{\mathrm{d}}t}}{{{\mathrm{d}}V}}\cdot \frac{{V\cos \gamma \sin \psi }}{{R\cos \phi }} $
$ \frac{{{\mathrm{d}}\phi }}{{{\mathrm{d}}V}} = \frac{{{\mathrm{d}}t}}{{{\mathrm{d}}V}}\cdot\frac{{V\cos \gamma \cos \psi }}{R} $
$ \begin{split}& \frac{{{\mathrm{d}}\psi }}{{{\mathrm{d}}V}} = \frac{{{\mathrm{d}}t}}{{{\mathrm{d}}V}}\left[ {\frac{{L\sin \sigma }}{{mV\cos \gamma }} + \frac{V}{R}\sin \psi \tan \phi \cos \gamma + 2{\omega _{\text{e}}}\sin \phi + } \right. \\&\qquad \left. {\frac{{R\omega _{\text{e}}^2\sin \psi \sin \phi \cos \phi }}{{V\cos \gamma }} - \frac{{2{\omega _{\text{e}}}\cos \psi \cos \phi \sin \gamma }}{{\cos \gamma }}} \right]\\[-1pt] \end{split} $
3.1　Longitudinal solutions
To analytically determine the altitude, FPA and bank angle, we define the altitude as a function of velocity.
$ H = f\left( V \right) = {a_6}{V^5} + {a_5}{V^4} + {a_4}{V^3} + {a_3}{V^2} + {a_2}V + {a_1} $
where ${a_i}(i = 1,2,3,4,5,6)$ are unknown coefficients.
The first-order and second-order derivatives of (18) can be expressed as
$ \frac{\text{d}H}{\text{d}V}={f}'\left(V\right)=5{a}_{6}{V}^{4}+4{a}_{5}{V}^{3}+3{a}_{4}{V}^{2}+2{a}_{3}V+{a}_{2}, $
$ \frac{{\text{d}}^{2}H}{\text{d}{V}^{2}}={f}''\left(V\right)=20{a}_{6}{V}^{3}+12{a}_{5}{V}^{2}+6{a}_{4}V+2{a}_{3}. $
By substituting (19) into (13), we can obtain the analytical expression of $\gamma $.
$ \gamma = \arcsin \left( {\frac{{{p_3}}}{{\sqrt {p_1^2 + p_2^2} }}} \right) - \arcsin \left( {\frac{{{p_2}}}{{\sqrt {p_1^2 + p_2^2} }}} \right) $
where
$ \left\{\begin{aligned}&{p}_{1}=V+\left(g-R{\omega }_{\text{e}}^{2}{\mathrm{cos}}^{2}\phi \right){f}'\left(V\right)\\&{p}_{2}=0.5R{\omega }_{\text{e}}^{2}\mathrm{sin}2\phi \mathrm{cos}\psi {f}'\left(V\right)\\&{p}_{3}=-\frac{{f}'\left(V\right)D}{m}\end{aligned}\right.. $
The first-order derivative of (21) can be expressed as
$\frac{\text{d}\gamma }{\text{d}V}=\frac{1}{{p}_{4}}\left[\frac{{f}'\left(V\right)}{m}\cdot\frac{\text{d}D}{\text{d}V}+\frac{D}{m}{f}''\left(V\right)+\\ \frac{\text{d}{p}_{1}}{\text{d}V}\mathrm{sin}\gamma +\frac{\text{d}{p}_{2}}{\text{d}V}\mathrm{cos}\gamma \right] $
where
$ \left\{\begin{aligned}&{p}_{4}={p}_{2}\mathrm{sin}\gamma -{p}_{1}\mathrm{cos}\gamma \\&\frac{\text{d}D}{\text{d}V}=\frac{\rho S}{2}\left(2V{C}_{D}+\frac{\partial {C}_{D}}{\partial Ma}\frac{{V}^{2}}{{V}_{\text{s}}}\right)\\&\frac{\text{d}{p}_{1}}{\text{d}V}=1-\left(\frac{2g}{R}+{\omega }_{\text{e}}^{2}{\mathrm{cos}}^{2}\phi \right){\left[{f}'\left(V\right)\right]}^{2}+R{\omega }_{\text{e}}^{\text{2}}\mathrm{sin}2\phi \cdot \\&\qquad {f}'\left(V\right)\frac{\text{d}\phi }{\text{d}V}+\left(g-R{\omega }_{\text{e}}^{\text{2}}{\mathrm{cos}}^{2}\phi \right){f}''\left(V\right)\\&\frac{\text{d}{p}_{2}}{\text{d}V}={\omega }_{\text{e}}^{\text{2}}\left\{\frac{\mathrm{sin}2\phi \mathrm{cos}\psi }{2}{\left[{f}'\left(V\right)\right]}^{2}+\left(\frac{R\mathrm{sin}2\phi \mathrm{cos}\psi }{2}\right)\cdot \right.\\&\qquad \left.{f}''\left(V\right)+R\left(\mathrm{cos}2\phi \mathrm{cos}\psi \frac{\text{d}\phi }{\text{d}V}-\frac{\mathrm{sin}2\phi \mathrm{sin}\psi }{2}\frac{\text{d}\psi }{\text{d}V}\right){f}'\left(V\right)\right\}\end{aligned} \right..$
Similarly, by substituting (23) into (14), the analytical expression of $\sigma $ can be expressed as
$ \sigma = \arcsin \left( {\frac{{{p_7} - {p_8}}}{{\sqrt {p_5^2 + p_6^2} }}} \right) - \arcsin \left( {\frac{{{p_6}}}{{\sqrt {p_5^2 + p_6^2} }}} \right) $
where
$ \left\{\begin{aligned}&{p}_{5}=\frac{LR{\omega }_{\text{e}}^{\text{2}}\mathrm{sin}\psi \mathrm{sin}2\phi {f}'\left(V\right)}{2mV{p}_{4}}\cdot\frac{\text{d}t}{\text{d}V}\\&{p}_{6}=\frac{L}{mV}\cdot\frac{\text{d}t}{\text{d}V}\\&{p}_{8}=\left[\left(\frac{V}{R}-\frac{g}{V}\right)\mathrm{cos}\gamma +2{\omega }_{\text{e}}\mathrm{cos}\phi \mathrm{sin}\psi +R{\omega }_{\text{e}}^{\text{2}}\left(\mathrm{cos}\gamma {\mathrm{cos}}^{2}\phi +\mathrm{cos}\psi \mathrm{sin}\gamma \frac{\mathrm{sin}2\phi }{2}\right)\right]\frac{\text{d}t}{\text{d}V}\\&{p}_{7}=\frac{1}{m{p}_{4}}\left(\frac{\text{d}D}{\text{d}V}{f}'\left(V\right)+D{f}''\left(V\right)\right)+\frac{\mathrm{sin}\gamma }{{p}_{4}}\left\{1-\left(\frac{2g}{R}+{\omega }_{\text{e}}^{\text{2}}{\mathrm{cos}}^{2}\phi \right){\left[{f}'\left(V\right)\right]}^{2}+\left(g-R{\omega }_{\text{e}}^{\text{2}}{\mathrm{cos}}^{2}\phi \right){f}''\left(V\right)+\right.\\&\qquad \left.\left(R{\omega }_{\text{e}}^{\text{2}}\mathrm{sin}2\phi \frac{\text{d}\phi }{\text{d}t}\frac{\text{d}t}{\text{d}V}\right){f}'\left(V\right)\right\}+\frac{{\omega }_{\text{e}}^{\text{2}}\mathrm{cos}\gamma }{{p}_{4}}\left\{\frac{R\mathrm{sin}2\phi \mathrm{cos}\psi }{2}{f}''\left(V\right)+\frac{\mathrm{sin}2\phi \mathrm{cos}\psi }{2}{\left[{f}'\left(V\right)\right]}^{2}+R\mathrm{cos}2\phi \mathrm{cos}\psi {f}'\left(V\right)\cdot \right.\\&\qquad \left.\frac{\text{d}\phi }{\text{d}t}\cdot\frac{\text{d}t}{\text{d}V}-\frac{R\mathrm{sin}2\phi \mathrm{sin}\psi }{2}\left(\frac{R{\omega }_{\text{e}}^{\text{2}}\mathrm{sin}2\phi \mathrm{sin}\psi }{2V\mathrm{cos}\gamma }+2{\omega }_{\text{e}}\mathrm{sin}\phi -2{\omega }_{\text{e}}\mathrm{cos}\phi \mathrm{cos}\psi \mathrm{tan}\gamma +\frac{V}{R}\mathrm{sin}\psi \mathrm{tan}\phi \mathrm{cos}\gamma \right){f}'\left(V\right)\frac{\text{d}t}{\text{d}V}\right\}\end{aligned}\right.. $
Finally, the altitude, FPA and bank angle are determined by the coefficients ${a_i}(i = 1,2,3,4,5,6)$.
3.2　Lateral solutions
Considering that (15)–(17) are highly nonlinear, it is necessary to simplify them. In this paper, we assume $\gamma \approx 0$ and $\dot \gamma \approx 0$, and ignore the $\omega _{\text{e}}^{\text{2}}$ term in (4) and (17). Meanwhile, since the HGV’s glide altitude ranges from 20 km to 80 km and $ {R_{\text{e}}} = 6371\;{\text{km}} $, it is evident that $H \ll {R_{\text{e}}}$. Therefore, during the glide phase, we assume that $R \approx {R^*} = {R_{\text{e}}} + ({H_0} + {H_{\text{f}}})/2$ and construct a nonlinear dynamic system $ {\boldsymbol{f}}({\boldsymbol{x}},V) $.
$ {\boldsymbol{f}}\left( {{\boldsymbol{x}},V} \right) = \left\{\begin{aligned}& {f_\lambda } = \frac{{{\mathrm{d}}\lambda }}{{{\mathrm{d}}V}} = - \frac{{mV\sin \psi }}{{D{R^*}\cos \phi }} \\& {f_\phi } = \frac{{{\mathrm{d}}\phi }}{{{\mathrm{d}}V}} = - \frac{{mV\cos \psi }}{{D{R^*}}} \\& {f_\psi } = \frac{{{\mathrm{d}}\psi }}{{{\mathrm{d}}V}} = - \frac{m}{{DV}}\left( {\frac{{L\sin \sigma }}{m} + 2V{\omega _{\text{e}}}\sin \phi } \right) - \\&\qquad \frac{{mV\sin \psi \tan \phi }}{{D{R^*}}} \end{aligned} \right. $
where ${H_0}$ represents the initial altitude, and $ {\boldsymbol{x}} = {\left[ {\lambda ,\phi ,\psi } \right]^{\text{T}}} $.
Then, we employ the Newton-Raphson method to solve the nonlinear dynamic system. The solutions of (27) obtained in the ith-iteration is denoted by $ {{\boldsymbol{x}}^i} = {[{\lambda ^i},{\phi ^i},{\psi ^i}]^{\text{T}}} $. As the Newton-Raphson method is an iterative method, there is
$ \boldsymbol{x}^{i+1}=\boldsymbol{x}^i+\Delta\boldsymbol{x}^i. $
By substituting (28) into (27), and then expanding (27) into a first-order Taylor series, we obtain
$ \boldsymbol{f}\left(\boldsymbol{x}^{i+1},V\right)=\boldsymbol{f}\left(\boldsymbol{x}^i,V\right)+\frac{\partial\boldsymbol{f}\left(\boldsymbol{x}^i,V\right)}{\partial\boldsymbol{x}^i}\left(\boldsymbol{x}^{i+1}-\boldsymbol{x}^i\right). $
Jacobian matrix is
$ \frac{{\partial {\boldsymbol{f}}\left( {{\boldsymbol{x}},V} \right)}}{{\partial {\boldsymbol{x}}}} = \left[ {\begin{array}{*{20}{c}} 0&{{{\partial {f_\lambda }} \mathord{\left/ {\vphantom {{\partial {f_\lambda }} {\partial \phi }}} \right. } {\partial \phi }}}&{{{\partial {f_\lambda }} \mathord{\left/ {\vphantom {{\partial {f_\lambda }} {\partial \psi }}} \right. } {\partial \psi }}} \\ 0&0&{{{\partial {f_\phi }} \mathord{\left/ {\vphantom {{\partial {f_\phi }} {\partial \psi }}} \right. } {\partial \psi }}} \\ 0&{{{\partial {f_\psi }} \mathord{\left/ {\vphantom {{\partial {f_\psi }} {\partial \phi }}} \right. } {\partial \phi }}}&{{{\partial {f_\psi }} \mathord{\left/ {\vphantom {{\partial {f_\psi }} {\partial \psi }}} \right. } {\partial \psi }}} \end{array}} \right] $
where
$ \left\{\begin{aligned} & \frac{\partial f_{\lambda}}{\partial\phi}=-\frac{mV\mathrm{tan}\phi\mathrm{sin}\psi}{DR^*\mathrm{cos}\phi}\text{，}\frac{\partial f_{\lambda}}{\partial\psi}=-\frac{mV\mathrm{cos}\psi}{DR^*\mathrm{cos}\phi} \\ & \frac{\partial f_{\phi}}{\partial\psi}=\frac{mV\mathrm{sin}\psi}{DR^*}\text{，}\frac{\partial f_{\psi}}{\partial\psi}=-\frac{mV\mathrm{cos}\psi\mathrm{tan}\phi}{DR^*} \\ & \frac{\partial f_{\psi}}{\partial\phi}=-\frac{2\omega_{\text{e}}m\mathrm{cos}\phi}{D}-\frac{mV\mathrm{sin}\psi\left(1+\mathrm{tan}^2\phi\right)}{DR^*}\end{aligned}\right.. $
Due to $ {{ - 2m{\omega _{\text{e}}}} \mathord{\left/ {\vphantom {{ - 2m{\omega _{\text{e}}}} D}} \right. } D} \leqslant 2.65 \times {10^{ - 4}} $ and ${{\partial {f_\psi }} \mathord{\left/ {\vphantom {{\partial {f_\psi }} {\partial \psi }}} \right. } {\partial \psi }} \leqslant 3 \times {10^{ - 4}}$, the first term $ {{ - 2m{\omega _{\text{e}}}\cos \phi } \mathord{\left/ {\vphantom {{ - 2m{\omega _{\text{e}}}\cos \phi } D}} \right. } D} $ in ${{\partial {f_\psi }} \mathord{\left/ {\vphantom {{\partial {f_\psi }} {\partial \phi }}} \right. } {\partial \phi }}$ and ${{\partial {f_\psi }} \mathord{\left/ {\vphantom {{\partial {f_\psi }} {\partial \psi }}} \right. } {\partial \psi }}$ can be ignored. Thus, from (29), there are
$\begin{split}& \frac{{{\text{d}}{\lambda ^{i + 1}}}}{{{\text{d}}V}} = d\left( V \right)\left[ {\frac{{\tan {\phi ^i}\sin {\psi ^i}}}{{\cos {\phi ^i}}}\left( {{\phi ^{i + 1}} - {\phi ^i}} \right) + } \right. \\&\qquad \left. {\frac{{\cos {\psi ^i}}}{{\cos {\phi ^i}}}\left( {{\psi ^{i + 1}} - {\psi ^i}} \right) + \frac{{\sin {\psi ^i}}}{{\cos {\phi ^i}}}} \right],\end{split} $
$ \frac{{{\text{d}}{\phi ^{i + 1}}}}{{{\text{d}}V}}\, = d\left( V \right)\left[ { - \sin {\psi ^i}\left( {{\psi ^{i + 1}} - {\psi ^i}} \right) + \cos {\psi ^i}} \right], $
$ \begin{split}& \frac{{{\text{d}}{\psi ^{i + 1}}}}{{{\text{d}}V}} = d\left( V \right)\left\{ {\left[ {\sin {\psi ^i}\left( {1 + {{\tan }^2}{\phi ^i}} \right)} \right]\left( {{\phi ^{i + 1}} - {\phi ^i}} \right) + } \right. \\& \quad\left. {\sin {\psi ^i}\tan {\phi ^i} + \frac{{{R^*}}}{{{V^2}}}\left( {\frac{{L\sin \sigma }}{m} + 2V{\omega _{\text{e}}}\sin {\phi ^i}} \right)} \right\}, \end{split} $
where $ d(V) = {{{ - mV} \mathord{\left/ {\vphantom {{ - mV} {DR}}} \right. } {(DR}}^*)} $.
Based on (33) and (34), we obtain
$ \left[ {\begin{array}{*{20}{c}} {\dfrac{{{\text{d}}{\phi ^{i + 1}}}}{{{\text{d}}V}}} \\ {\dfrac{{{\text{d}}{\psi ^{i + 1}}}}{{{\text{d}}V}}} \end{array}} \right] = d\left( V \right)\left\{ {{\boldsymbol{A}}\left[ {\begin{array}{*{20}{c}} {{\phi ^{i + 1}}} \\ {{\psi ^{i + 1}}} \end{array}} \right] + \left( {{{\boldsymbol{B}}_1} + {{\boldsymbol{B}}_2}} \right)} \right\} $
where
$ \left\{ \begin{aligned}& {\boldsymbol{A}} = \left[ {\begin{array}{*{20}{c}} 0&{ - \sin {\psi ^i}} \\ {\sin {\psi ^i}\left( {1 + {{\tan }^2}{\phi ^i}} \right)}&0 \end{array}} \right]\quad \\& {{\boldsymbol{B}}_1} = \dfrac{{{R^*}}}{{{V^2}}}\left[ {\begin{array}{*{20}{c}} 0 \\ {\left( {\dfrac{{L\sin \sigma }}{m} + 2V{\omega _{\text{e}}}\sin {\phi ^i}} \right)} \end{array}} \right] \\& {{\boldsymbol{B}}_2} = \left[ {\begin{array}{*{20}{c}} {\sin {\psi ^i}{\psi ^i} + \cos {\psi ^i}} \\ { - \sin {\psi ^i}\left( {1 + {{\tan }^2}{\phi ^i}} \right){\phi ^i} + \sin {\psi ^i}\tan {\phi ^i}} \end{array}} \right]\end{aligned} \right.. $
Multiplying both sides of (35) by $ {\boldsymbol{\varTheta}} = \exp \left( \displaystyle\int_{{V_0}}^V - {\boldsymbol{A}}\cdot d\left( u \right){\text{d}}u \right) $ yields
$ {\boldsymbol{\varTheta}} \left[ {\begin{array}{*{20}{c}} {\dfrac{{{\text{d}}{\phi ^{i + 1}}}}{{{\text{d}}V}}} \\ {\dfrac{{{\text{d}}{\psi ^{i + 1}}}}{{{\text{d}}V}}} \end{array}} \right] - {\boldsymbol{\varTheta}} \cdot d\left( V \right) \cdot {\boldsymbol{A}}\left[ {\begin{array}{*{20}{c}} {{\phi ^{i + 1}}} \\ {{\psi ^{i + 1}}} \end{array}} \right] = {\boldsymbol{\varTheta}} \cdot d\left( V \right)\left[ {{{\boldsymbol{B}}_1} + {{\boldsymbol{B}}_2}} \right] $
where ${V_0}$ is the initial velocity.
By integrating (37), we obtain the analytical expressions for $\phi $ and $\psi $.
$ \begin{split}& \left[ {\begin{array}{*{20}{c}} {{\phi ^{i + 1}}} \\ {{\psi ^{i + 1}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\phi _0}} \\ {{\psi _0}} \end{array}} \right] + \int_{{V_0}}^V {d\left( u \right){\boldsymbol{\varTheta}} \cdot {{\boldsymbol{B}}_1}{\text{d}}u} + \left. {\left[ { - {\boldsymbol{\varTheta}} \cdot {{\boldsymbol{A}}^{ - 1}}} \right]} \right|_{{V_0}}^V \cdot \\&\qquad {{\boldsymbol{B}}_2}= {{\boldsymbol{\varTheta}} ^{ - 1}}\left[ {\begin{array}{*{20}{c}} {{\phi _0}} \\ {{\psi _0}} \end{array}} \right] + {{\boldsymbol{\varTheta}} ^{ - 1}} \cdot \int_{{V_0}}^V {d\left( u \right){\boldsymbol{\varTheta}} \cdot {{\boldsymbol{B}}_1}{\text{d}}u} + \\&\qquad\qquad\qquad \left[ { - {\boldsymbol{\varTheta}} + {{\boldsymbol{I}}_{2 \times 2}}} \right] \cdot {{\boldsymbol{A}}^{ - 1}} \cdot {{\boldsymbol{B}}_2}\\[-1pt] \end{split} $
where $ {{\boldsymbol{I}}_{2 \times 2}} $ is the identity matrix, ${\phi _0}$ is the initial latitude, and ${\psi _0}$ is the initial heading angle.
By using the obtained solutions for $H$, $\gamma $, $\phi $, and $\psi $, we can reduce the order of the reentry motion equations. In these solutions, $V$ serves as the independent variable. Therefore, the reentry motion equations can be simplified to first order, involving only ${{{\text{d}}\lambda } \mathord{\left/ {\vphantom {{{\text{d}}\lambda } {{\text{d}}V}}} \right. } {{\text{d}}V}}$.
4.　Online guidance algorithm
Less
4.1　Calculation of AOA
We construct the optimal AOA profile as a piecewise linear function of velocity, which is specifically expressed as follows:
$ \alpha \left( V \right) = \left\{\begin{aligned}& {\alpha _{{\text{max}}}},\; V \gt {V_{{\alpha _1}}} \\& \frac{{{\alpha _{L/{D_{\max }}}} - {\alpha _{{\text{max}}}}}}{{{V_{{\alpha _2}}} - {V_{{\alpha _1}}}}}\left( {V - {V_{{\alpha _1}}}} \right) + {\alpha _{{\text{max}}}},\;{V_{{\alpha _2}}} \leqslant V \leqslant {V_{{\alpha _1}}} \\& {\alpha _{L/{D_{\max }}}},\;V \lt {V_{{\alpha _2}}} \end{aligned} \right. $
where $ {\alpha _{L/{D_{\max }}}} $ is the maximum lift-to-drag ratio AOA, $ {V_{{\alpha _1}}} $ and $ {V_{{\alpha _2}}} $ are two given velocity nodes, respectively.
4.2　Calculation of altitude-velocity profile
Since ${H_0}$, ${H_{\text{f}}}$, ${V_0}$ and ${V_{\text{f}}}$ are known, two equations can be obtained by substituting them into (18).
$ {H_0} = {a_6}V_0^5 + {a_5}V_0^4 + {a_4}V_0^3 + {a_3}V_0^2 + {a_2}{V_0} + {a_1} $
$ {H_{\text{f}}} = {a_6}V_{\text{f}}^5 + {a_5}V_{\text{f}}^4 + {a_4}V_{\text{f}}^3 + {a_3}V_{\text{f}}^2 + {a_2}{V_{\text{f}}} + {a_1} $
Similarly, ${\gamma _0}$ and ${\gamma _{\text{f}}}$ are also known. Assuming that the Earth’s rotation term in (4) is ignored, we can derive two equations according to (13) and (19).
$ \frac{{m{V_0}\sin {\gamma _0}}}{{ - {D_0} - mg\sin {\gamma _0}}} = 5{a_6}V_0^4 + 4{a_5}V_0^3 + 3{a_4}V_0^2 + 2{a_3}{V_0} + {a_2} $
$ \frac{{m{V_{\text{f}}}\sin {\gamma _{\text{f}}}}}{{ - {D_{\text{f}}} - mg\sin {\gamma _{\text{f}}}}} = 5{a_6}V_{\text{f}}^4 + 4{a_5}V_{\text{f}}^3 + 3{a_4}V_{\text{f}}^2 + 2{a_3}{V_{\text{f}}} + {a_2} $
To accurately determine the value of ${a_i}(i = 1,2,3, 4,5,6)$, we require two additional equations. Considering that the velocity of HGV decreases monotonically during the glide phase, we select two velocity nodes, ${V_1} = {V_0} - {{\left( {{V_0} - {V_{\text{f}}}} \right)} \mathord{\left/ {\vphantom {{\left( {{V_0} - {V_{\text{f}}}} \right)} 3}} \right. } 3}$ and ${V_2} = {V_0} - {{2\left( {{V_0} - {V_{\text{f}}}} \right)} \mathord{\left/ {\vphantom {{2\left( {{V_0} - {V_{\text{f}}}} \right)} 3}} \right. } 3}$, for analysis. The altitudes corresponding to ${V_1}$ and ${V_2}$ are represented by ${H_1}$ and ${H_2}$, respectively. They can be optimized based on the IPSO algorithm, which is described in detail in Subsection 4.3.
Thus, there are
$ {H_1} = {a_6}V_1^5 + {a_5}V_1^4 + {a_4}V_1^3 + {a_3}V_1^2 + {a_2}{V_1} + {a_1}, $
$ {H_2} = {a_6}V_2^5 + {a_5}V_2^4 + {a_4}V_2^3 + {a_3}V_2^2 + {a_2}{V_2} + {a_1}. $
From (40)−(45), there are
$ {\boldsymbol{U}} = {{\boldsymbol{Y}}^{ - 1}}{\boldsymbol{Z}} $
where
$\left\{\begin{aligned}&{\boldsymbol{U}} = \left[ {\begin{array}{*{20}{c}} {{a_6}} \\ {{a_5}} \\ {{a_4}} \\ {{a_3}} \\ {{a_2}} \\ {{a_1}} \end{array}} \right]\\&{\boldsymbol{Y}} = \left[ {\begin{array}{*{20}{c}} {V_0^5}&{V_0^4}&{V_0^3}&{V_0^2}&{{V_0}}&1 \\ {V_1^5}&{V_1^4}&{V_1^3}&{V_1^2}&{{V_1}}&1 \\ {V_2^5}&{V_2^4}&{V_2^3}&{V_2^2}&{{V_2}}&1 \\ {V_{\text{f}}^5}&{V_{\text{f}}^4}&{V_{\text{f}}^3}&{V_{\text{f}}^2}&{{V_{\text{f}}}}&1 \\ {5V_0^4}&{4V_0^3}&{3V_0^2}&{2{V_0}}&1&0 \\ {5V_{\text{f}}^4}&{4V_{\text{f}}^3}&{3V_{\text{f}}^2}&{2{V_{\text{f}}}}&1&0 \end{array}} \right]\\&{\boldsymbol{Z}} = {\left[ {{H_0},{H_1},{H_2},{H_{\text{f}}},\frac{{m{V_0}\sin {\gamma _0}}}{{ - {D_0} - mg\sin {\gamma _0}}},\frac{{m{V_{\text{f}}}\sin {\gamma _{\text{f}}}}}{{ - {D_{\text{f}}} - mg\sin {\gamma _{\text{f}}}}}} \right]^{\text{T}}}\end{aligned}\right.. $
4.3　IPSO
In PSO, every particle represents a possible solution and is associated with a position vector and a velocity vector. Assume that the number of particles in the swarm is Np, and the number of optimization parameters is Nop. The position vector of the Ith particle is denoted by ${{\boldsymbol{X}}_I} = ({x_{I,1}},{x_{I,2}}, \cdots ,{x_{I,{N_{{\text{op}}}}}})$, and the velocity vector is denoted by $ {{\boldsymbol{V}}_I} = ({v_{I,1}},{v_{I,2}}, \cdots ,{v_{I,{N_{{\text{op}}}}}}) $. The Ith best particle is denoted by $ {{\boldsymbol{P}}_I} = ({p_{I,1}},{p_{I,2}}, \cdots ,{p_{I,{N_{{\text{op}}}}}}) $, and the best particle in the swarm is denoted by $ {{\boldsymbol{P}}_{\text{b}}} = ({p_{{\text{b}},1}},{p_{{\text{b}},2}}, \cdots ,{p_{{\text{b}},{N_{{\text{op}}}}}}) $. The swarm evolves as follows [24]:
$ v_{I,j}^{t_{\text{ci}}+1}=wv_{I,j}^{t_{\text{ci}}}+c_1r_1\left(p_{I,j}-x_{I,j}^{t_{\text{ci}}}\right)+c_2r_2\left(p_{\text{b},j}-x_{I,j}^{t_{\text{ci}}}\right), $
$ x_{I,j}^{t_{\text{ci}}+1}=x_{I,j}^{t_{\text{ci}}}+v_{I,j}^{t_{\text{ci}}+1}, $
where ${t_{{\text{ci}}}}$ is the current iteration number, $w$ is the inertia weight. ${c_1}$ and ${c_2}$ are cognitive factors, in general, $ {c_1} = {c_2} = 2.0 $. $ {r_1} $ and $ {r_2} $ are independent random numbers between 0 and 1.
Since the PSO algorithm cannot solve constrained problems and suffers from local convergence issues, several improvements have been proposed to enhance the performance of the PSO algorithm. The improvements are detailed later.
4.3.1　Objective function
We utilize the penalty function method to address the reentry constraints. Therefore, we summarize the objective function as follows:
$ \min J = \left| {{R_{{\text{stogo}}}}\left( {{t_{\text{f}}}} \right) - {R_{{\text{stogof}}}}} \right| + {W_1}\left( {{P_1} + {P_2} + {P_3} + {P_4}} \right) $
where $ W_1=1\ 000 $ is the penalty parameter.
$ \left\{ \begin{gathered} {P_1} = \max (0,\dot Q - {{\dot Q}_{\max }}) \\ {P_2} = \max (0,n - {n_{\max }}) \\ {P_3} = \max (0,q - {q_{\max }}) \\ {P_4} = \max (0,\left| \sigma \right| - {\sigma _{\max }}) \\\end{gathered} \right. $
4.3.2　Restriction of optimization parameter boundaries
In practical problems, the optimization parameters are always restricted by numerical boundaries. Therefore, the position and the velocity of a particle should be constrained.
In PSO, expanding the parameter optimization range increases the search space, reduces convergence speed, raises computational costs, and makes parameter adjustment more difficult. For this study, generally, $ x_{{\text{up}}}^j = {H_0} $, and $ x_{{\text{low}}}^j = {H_{\text{f}}} $. To ensure that the reentry trajectory adheres to path constraints while also balancing search efficiency and computational resource consumption, it is necessary to further narrow the range of optimization parameters.
$ x_{{\text{low}}}^j = \max \left( {{H_{\dot Q}},{H_n},{H_q}} \right) $
where, according to (8), $ ({H_{\dot Q}},{H_n},{H_q}) $ can be expressed as
$ \left\{ \begin{gathered} {H_{\dot Q}} \geqslant \frac{2}{\beta }\ln \frac{{{k_{\dot Q}}\rho _0^{0.5}{V^{3.15}}}}{{{{\dot Q}_{\max }}}} \\ {H_n} \geqslant \frac{1}{\beta }\ln \frac{{{\rho _0}{V^2}}}{{2{q_{\max }}}} \\ {H_q} \geqslant \frac{1}{\beta }\ln \frac{{\sqrt {C_D^2 + C_L^2} {\rho _0}{V^2}S}}{{2{n_{\max }}mg}} \\ \end{gathered} \right.. $
Then, to address the issue of particles frequently moving around the optimal position without precisely converging as they approach the optimal solution, we adopt a strategy to reduce their velocity. Let the velocity boundary shrink linearly with iteration time.
$ \left\{ \begin{gathered} v_{{\text{up}}}^j\left( {{t_{{\text{ci}}}}} \right) = \left( {x_{{\text{up}}}^j - x_{{\text{low}}}^j} \right) - \frac{{\left( {x_{{\text{up}}}^j - x_{{\text{low}}}^j} \right)}}{{t_{{\text{ci}}}^{{\text{max}}}}}\left( {{t_{{\text{ci}}}} - 1} \right) \\ v_{{\text{low}}}^j\left( {{t_{{\text{ci}}}}} \right) = v_{{\text{up}}}^j\left( {{t_{{\text{ci}}}}} \right) \\ \end{gathered} \right. $
where $ t_{{\text{ci}}}^{{\text{max}}} $ is the maximum iteration number.
4.3.3　Inertia weight and cognitive factors adjustment strategy
The inertia weight reflects the effect of the velocity on a particle’s movement and influences the evolution. A large weight helps global searching, while a small weight benefits convergence. Here, the weight is assigned with a large value initially in order to fully detect the searching space; as the swarm evolves, the weight decreases to facilitate convergence. Here, a linear decreasing inertia weight is utilized.
$ w = \left( {{w_1} - {w_0}} \right)\frac{{{t_{{\text{ci}}}} - 1}}{{t_{{\text{ci}}}^{{\text{max}}}}} + {w_0} $
where $0 \lt w \lt 1$, $ {w_0} $ is the initial value of the weight and $ {w_1} $ is the final value.
$ {c_1} $ enhances the local search ability of particles, making them more inclined to explore their own historical optimal solutions, while $ {c_2} $ promotes the movement of particles towards the swarm’s optimal solution and enhances global search capabilities. Reasonably adjusting these parameters can help achieve a balance between local and global searches in the algorithm, thereby improving the overall optimization effect.
In addition, if $ {c_1} $ is too high, it may cause the algorithm to prematurely converge to a local optimum. Conversely, if $ {c_2} $ is too high, it will increase the randomness of the search process, thus affecting the convergence speed. Therefore, a linear adjustment strategy is adopted.
$ \left\{ \begin{gathered} {c_1} = c_1^0 - 1.1\frac{{{t_{{\text{ci}}}} - 1}}{{t_{{\text{ci}}}^{{\text{max}}}}} \\ {c_2} = c_2^0 + \frac{{{t_{{\text{ci}}}} - 1}}{{t_{{\text{ci}}}^{{\text{max}}}}} \\ \end{gathered} \right. $
where $ c_1^0 $ and $ c_2^0 $ are the initial values of $ {c_1} $ and $ {c_2} $, respectively.
4.4　Bank reversal strategy
For lateral guidance, the heading angle error corridor is used to determine the bank angle symbol. The heading angle error corridor is
$ \left| {\Delta {\psi _{{\text{th}}}}} \right| = \left\{\begin{aligned}&10,\; 6\;000 \lt V \leqslant {V_0} \\&15,\; 3\;000 \lt V \leqslant 6\;000 \\& \frac{7}{{1\;200}}\left( {V - 3\;000} \right) + 15,\; 1\;800 \lt V \leqslant 3\;000 \\&8,\; V \leqslant {V_{\text{f}}}\end{aligned} \right.. $
Once the difference $\Delta \psi $ between the heading angle and the line of sight angle exceeds the heading angle error corridor boundary, the bank angle symbol is changed. The symbol transformation strategy is
$ {\text{sgn}}\left( \sigma \right) = \left\{\begin{aligned}& - 1,\;\;\Delta \psi \geqslant \left| {\Delta {\psi _{{\text{th}}}}} \right| \\& {\text{sgn}}\left( {{\sigma _{\text{p}}}} \right),\;\; - \left| {\Delta {\psi _{{\text{th}}}}} \right| \lt \Delta \psi \lt \left| {\Delta {\psi _{{\text{th}}}}} \right| \\&1,\;\;\Delta \psi \leqslant - \left| {\Delta {\psi _{{\text{th}}}}} \right| \end{aligned}\right. $
where $ {\sigma _{\text{p}}} $ is the bank angle symbol of the previous step, $ \Delta \psi = \psi - {\psi _{{\text{los}}}} $, $ {\psi _{{\text{los}}}} $ is expressed as follows:
$ \tan {\psi _{{\text{los}}}} = \frac{{\sin \left( {{\lambda _{\text{T}}} - \lambda } \right)}}{{\cos \phi \tan {\phi _{\text{T}}} - \sin \phi \cos \left( {{\lambda _{\text{T}}} - \lambda } \right)}}. $
4.5　Algorithm framework
Fig. 1 shows the framework diagram of the online reentry guidance algorithm. A brief description of the guidance process is as follows:
Step 1 Set the initial states and optimize the altitude profile (see Step 3).
Step 2 Calculate the AOA with (39).
Step 3 Update the altitude profile.
Step 3.1 Set the IPSO parameters and ${t_{{\text{ci}}}} = 0$.
Step 3.2 Initialize the swarm and get ${{\boldsymbol{P}}_I}$ and ${{\boldsymbol{P}}_{\text{b}}}$.
Step 3.3 Update $w$, ${c_1}$ and ${c_2}$ with (56) and (57).
Step 3.4 Update the particle position and velocity, and then update ${{\boldsymbol{P}}_I}$ and ${{\boldsymbol{P}}_{\text{b}}}$.
Step 3.5 Update the velocity boundary with (55), and then restrict the particle position and velocity.
Step 3.6 Calculate ${R_{{\text{stogo}}}}$. If ${R_{{\text{stogo}}}} \geqslant {R_{{\text{stogof}}}}$, set ${t_{{\text{ci}}}} = {t_{{\text{ci}}}} + 1$ and return to Step 3.3. Otherwise, output ${H_1}$ and ${H_2}$, and calculate ${a_i}(i = 1,2,3,4,5,6)$ with (46).
Step 4 Update the flight states.
Step 4.1 Calculate the altitude, FPA and bank angle with (18), (21) and (25).
Step 4.2 Calculate the latitude and heading angle with (38), and then integral ${\mathrm{d}}\lambda/{\mathrm{d}}V $ in (17) to calculate the longitude.
Step 5 If $V \geqslant {V_{\text{f}}}$, set $V = V - 5$ and return to Step 2. Otherwise, end the algorithm.
5.　Simulation tests
Less
In this paper, we conduct a series of simulation tests using the CAV-H [25]—a reentry vehicle concept model developed by Lockheed Martin. All simulations are executed on a personal computer equipped with a 2.2 GHz CPU and 16 GB of RAM. The reference area of the CAV-H is 0.484 m2, and the mass is 907.2 kg.
5.1　HGT analytical solutions performance analysis
To evaluate the performance of the HGT analytical solutions derived in this paper regarding computational accuracy and computational complexity, this subsection compares them with the original reentry motion equations (ORMEs) and the reduced-order motion equations (ROMEs) proposed in [20]. Specifically, trajectories are generated based on these motion equations, and their computational accuracy and efficiency are compared. For the differential equations involved in ROMEs and ORMEs, this subsection employs the Runge-Kutta (RK) method for numerical integration.
This subsection considers three examples, with the initial positions of the vehicle all set to ${\lambda _0} = - 72^\circ $, ${\phi _0} = - 20^\circ $. The three examples utilize different orders of altitude-velocity profiles, specifically the 5th order, the 4th order, and the 2nd order, with their specific forms illustrated in Fig. 2.
Fig. 3 illustrates the ground trajectories generated based on the HGT analytical solutions, ROMEs, and ORMEs across all examples. It is evident that the trajectory generated using the HGT analytical solutions is the closest to the result obtained from ORMEs, with significantly higher computational accuracy compared to ROMEs. This improvement in accuracy is attributed to the comprehensive consideration of the Earth’s rotation during the derivation of the HGT analytical solutions.
Table 1 presents the computation times required to generate trajectories using the three motion equations. The results indicate that the computation time for the trajectory generated by the HGT analytical solutions is within 10 ms, demonstrating a significantly higher computational efficiency compared to ORMEs and ROMEs. By reducing the order of the reentry motion equations to the first order, the computational complexity is greatly diminished, thereby significantly enhancing computational efficiency.
5.2　Online guidance algorithm performance analysis
In this subsection, we select four examples with different targets to verify and analyze the proposed algorithm. Table 2 gives the simulation conditions for these examples, while the parameter settings for the IPSO algorithm are detailed in Table 3. The path constraints and control constraints are set to: ${\dot Q_{\max }} = 5\;000\;{\text{kW/}}{{\text{m}}^2}$, ${n_{\max }} = 5g$, ${q_{\max }} = 120\;{\text{kPa}}$, $ {\alpha _{\max }} = 25^\circ $, and $ {\sigma _{\max }} = 80^\circ $.
5.2.1　Different examples simulations and discussion
This subsection verifies the effectiveness and adaptability of the proposed algorithm across different examples. The simulation results are shown in Fig. 4-Fig. 11 and Table 4. The results demonstrate that the reentry constraints are consistently well satisfied in each example. Because the designed altitude-velocity profile considers the terminal altitude, velocity, and FPA requirements, the simulation results of the proposed algorithm satisfy the terminal constraints very well. However, an error of 0.01° in the terminal FPA is present only in Example 4. This error is negligible and may be attributable to the simplification of the reentry motion equations during the derivation of the HGT analytical solutions. In addition, although the PSO algorithm faces challenges when dealing with constraints, the IPSO algorithm has proven effective in addressing these issues. It successfully determines the optimal altitude-velocity profile that satisfies all constraints, ensuring that the range-to-go of the vehicle is within 200 km.
5.2.2　Comparison of different algorithms
This subsection utilizes the simulation conditions described in Example 1 to explore three different guidance algorithms. The first is the proposed algorithm in this paper. The second involves combining the basic PSO algorithm (BPSO) and the RK method into the presented online guidance algorithm framework, and this algorithm is called BPSO+RK. Finally, the PSO algorithm proposed by Zhou [21], along with the RK method, is integrated into the online guidance algorithm framework. This algorithm is designated as ZPSO+RK.
Fig. 12 and Table 5 show the simulation results for three algorithms used in optimizing the altitude-velocity profile. The results indicate that the IPSO algorithm exhibits the fastest convergence speed, achieving convergence after just 18 iterations and completing the optimization process in only 0.87 s. Given that the proposed algorithm employs the HGT analytical solutions when predicting the range-to-go, it significantly reduces the computational burden compared to the integration of reentry motion equations using the RK method. Coupled with the IPSO algorithm’s rapid convergence, the time required to optimize a single profile is controlled to be within 2 s, meeting the requirements for online guidance. Although the BPSO achieves the highest convergence accuracy, it converges the slowest, requiring 44 iterations and taking 5.3 s for the optimization process. In comparison, the ZPSO improves convergence speed to some extent, limiting the number of iterations to 30, but the time required to optimize the profile still reaches 4.1 s. Moreover, the ZPSO sacrifices some convergence accuracy, indicating a need for further enhancement.
Through the analysis of Fig. 13, it has been observed that the three algorithms can guide the vehicle to within 200 km of the target while satisfying the reentry constraints. Among these algorithms, the range-to-go for the proposed algorithm is 21.07 km, for the BPSO+RK it is 15.95 km, and for the ZPSO+RK it is 103.08 km. Although the range-to-go of the proposed algorithm is not the smallest, it satisfies the terminal range requirements. The computation times of the three algorithms are shown in Table 6. The computation time of the proposed algorithm is about 8.3 s. Considering computational efficiency, the proposed algorithm is more suitable for the online guidance than the other two algorithms.
5.3　Monte Carlo simulations
To further evaluate the computational accuracy, efficiency, and robustness of the proposed algorithm, this subsection examines the perturbations of the aerodynamic lift and drag coefficients, with the perturbation amplitude constrained within ±10%. A total of 100 groups of perturbation conditions are randomly selected for the Monte Carlo simulations. The initial and terminal conditions of the simulation are aligned with those specified in Example 1 in Table 2.
The simulation results are shown in Fig. 14, where Fig. 14(a) and Fig. 14(b) illustrates the flight trajectories and the ground trajectories, respectively. It is observed that under the disturbance conditions, the trajectories remain relatively smooth. Fig. 14(c) illustrates the altitude-velocity profiles under these disturbance conditions. The results indicate that the trajectories effectively satisfy the quasi-equilibrium glide condition and are evenly distributed within the reentry corridor, thereby ensuring compliance with the path constraints. Fig. 14(d) presents the FPA-velocity profiles under the disturbance conditions, demonstrating that all FPAs satisfy the terminal constraint conditions. Fig. 14(e) shows the bank angle-velocity profiles under the same conditions, with results indicating that all bank angles adhere to the control constraints. Fig. 14(f) illustrates the distribution of the terminal landing point, revealing that the range-to-go is effectively controlled within 200 km, thus meeting the mid-terminal handover conditions. Finally, Fig. 14(g) illustrates the computation time distribution under the disturbance conditions, with the computation time is controlled in 7−10 s, further validating the high computational efficiency of the proposed algorithm. In summary, the simulation results comprehensively demonstrate that the proposed algorithm exhibits high computational accuracy, efficiency, and robustness under disturbance conditions.
6.　Conclusions
Less
This paper proposes an online guidance algorithm based on the HGT analytical solutions to address the challenge of integral-calculated trajectories that often fail to satisfy the requirements for online guidance. Specifically, we derive the analytical solutions for the longitudinal trajectory and the recursive analytical solutions for the lateral trajectory. These solutions effectively reduce the order of the reentry motion equations to the first order. On this basis, we introduce a profile updating strategy that replaces traditional trajectory tracker design. This strategy employs a fifth-order polynomial profile, which inherently satisfies the constraints of terminal altitude, velocity, and FPA. Furthermore, we improve the PSO algorithm by narrowing the parameter optimization range and adaptively adjusting the inertia weight and cognitive factors. This significantly accelerates the algorithm’s global convergence. By combining the profile updating strategy with the IPSO algorithm, we present an online guidance algorithm capable of achieving effective vehicle guidance. This is accomplished by optimizing the altitude-velocity profile for each guidance cycle and incorporating a bank reversal strategy. The proposed algorithm does not require multiple iterations or extensive integral computations for online guidance. The time required to update guidance commands within a single guidance cycle is less than 2 s, showcasing a significant improvement in computational efficiency compared to traditional algorithms.
It should be pointed out that the analytical solutions derived in this paper have not been fully resolved and still involve integral terms, which to some extent reduces the speed and effectiveness of the online guidance. In future work, we could consider using the Chebyshev series or Euler series to approximate the integral terms involved in the analytical solutions with high precision. This would help achieve a fully analytical calculation of HGT and further enhance the online capability of the reentry guidance algorithm. Additionally, exploring more advanced optimization algorithms, such as neural networks and deep learning, for the design of longitudinal and lateral guidance strategies can also further improve guidance efficiency.
References
Less
1
DING Y B, YUE X K, CHEN G S, et al. Review of control and guidance technology on hypersonic vehicle. Chinese Journal of Aeronautics, 2022, 35(7): 1–18.
2
ZHANG Y L, XIE Y. Review of trajectory planning and guidance methods for gliding vehicles. Acta Aeronautica et Astronautica Sinica, 2020, 41(1): 50–62. (in Chinese)
3
GUO J, ZHENG J K, WANG H N, et al. Review of research on reentry guidance methods and hot issues of hypersonic gliding vehicle. Aerospace Technology, 2022, 2022(1): 54–63.
4
TIAN B L, LI Z Y, WU S Y, et al. Reentry trajectory optimization, guidance and control methods for reusable launch vehicles: review. Acta Aeronautica et Astronautica Sinica, 2020, 41(11): 624072. (in Chinese)
5
HE R Z, LIU L H, TANG G J, et al. Entry trajectory generation without reversal of bank angle. Aerospace Science and Technology, 2017, 71: 627–635.
6
MORIO V, CAZAURANG F, VERNIS P. Flatness-based hypersonic reentry guidance of a lifting-body vehicle. Control Engineering Practice, 2009, 17(5): 588–596.
7
YAN H, HE Y Z. Drag-tracking guidance for entry vehicles without drag rate measurement. Aerospace Science and Technology, 2015, 43: 372–380.
8
MEASE K D, CHEN D T, TEUFEL P, et al. Reduced-order entry trajectory planning for acceleration guidance. Journal of Guidance, Control, and Dynamics, 2002, 25(2): 257–266.
9
LEAVITT J A, MWASE K D. Feasible trajectory generation for atmospheric entry guidance. Journal of Guidance, Control, and Dynamics, 2007, 30(2): 473–481.
10
WANG X, TANG S J, QI S. Predictor-corrector entry guidance with terminal altitude constraint. Tactical Missile Technology, 2018, 4: 70–77. (in Chinese)
11
WANG T, ZHANG H B, TANG G J. Predictor-corrector entry guidance with waypoint and no-fly zone constraints. Acta Astronautica, 2017, 138: 10–18.
12
WANG T, ZHANG H B, ZENG L, et al. A robust predictor-corrector entry guidance. Aerospace Science and Technology, 2017, 66: 103–111.
13
ZHANG P, DU Y L, XIANG K. Constrained predictive-corrector reentry guidance for high lift-to-drag RLV. Flight Dynamics, 2018, 36(3): 70–74. (in Chinese)
14
XU H, CAI G B, MU C X, et al. Analytical reentry guidance framework based on swarm intelligence optimization and altitude-energy profile. Chinese Journal of Aeronautics, 2023, 36(12): 336–348.
15
PAN L, PENG S C, XIE Y, et al. 3D guidance for hypersonic reentry gliders based on analytical prediction. Acta Astronautica, 2020, 167: 42–51.
16
YU W B, CHEN W C. Entry guidance with real-time planning of reference based on analytical solutions. Advances in Space Research, 2015, 55(9): 2325–2345.
17
YU W B, CHEN W C, JIANG Z G, et al. Analytical entry guidance based on pseudo-aerodynamic profiles. Aerospace Science and Technology, 2017, 66: 315–331.
18
YU W B, YANG J, CHEN W C. Entry guidance based on analytical trajectory solutions. IEEE Trans. on Aerospace and Electronic Systems, 2021, 58(3): 2438–2466.
19
YANG J, YU W B, CHEN W C, et al. Chebyshev-series solutions for nonlinear systems with hypersonic gliding trajectory example. Aerospace Science and Technology, 2023, 140: 108424.
20
ZENG L, ZHANG H B, ZHENG W. A three-dimensional predictor-corrector entry guidance based on reduced-order motion equations. Aerospace Science and Technology, 2018, 73: 223–231.
21
ZHOU H Y, WANG X G, CUI N G. A novel reentry trajectory generation method using improved particle swarm optimization. IEEE Trans. on Vehicular Technology, 2019, 68(4): 3212–3223.
22
ZHOU H Y, WANG X G, CUI N G. Glide guidance for reusable launch vehicles using analytical dynamics. Aerospace Science and Technology, 2020, 98: 105678.
23
ZHOU H Y, LI X, BAI Y L, et al. Optimal guidance for hypersonic vehicle using analytical solutions and an intelligent reversal strategy. Aerospace Science and Technology, 2023, 132: 108053.
24
SOO K K, SIU Y M, CHAN W S, et al. Particle-swarm-optimization-based multiuser detector for CDMA communications. IEEE Trans. on Vehicular Technology, 2007, 56(5): 3006–3013.
25
PHILLIPS T. A common aero vehicle (CAV) model, description, and employment guide. Schafer Corporation for AFRL and AFSPC, Arlington, VA 2003.
Appendix
Less
Year 2026 volume 37 Issue 3
PDF
79
43
Cite this Article
BibTeX
Article Info
doi: 10.23919/JSEE.2026.000089
  • Receive Date：2024-08-26
  • Online Date：2026-08-14
  • Published：2026-06-18
Article Data
Affiliations
History
  • Received：2024-08-26
Affiliations
    1Control and Simulation Center, Harbin Institute of Technology, Harbin 150001, China
    2National Key Laboratory of Modeling and Simulation for Complex Systems, Harbin 150001, China

Corresponding:

CHAO Tao
References
Share
https://castjournals.cast.org.cn/joweb/jsee/EN/10.23919/JSEE.2026.000089
Share to
QR

Scan QR to access full text

Cite this article
BibTeX
Citations
表12种不同金属材料的力学参数
科
Family		 属数
Number of
genus		 种数
Number of
species		 占总种数比例
Percentage of
total species (%)		 属
Genus		 种数
Number of
species		 占总种数比例
Percentage of total
species (%)
鹅膏菌科Amanitaceae 2 11 5.26 鹅膏菌属 Amanita 10 4.78
小菇科 Mycenaceae 2 12 5.74 丝盖伞属 Inocybe 5 2.39
多孔菌科 Polyporaceae 8 14 6.70 蜡蘑属 Laccaria 5 2.39
红菇科 Russulaceae 3 23 11.00 小皮伞属 Marasmius 6 2.87
小菇属 Mycena 11 5.26
光柄菇属 Pluteus 5 2.39
红菇属 Russula 17 8.13
栓菌属 Trametes 5 2.39
关闭全屏
  • BibTeX
  • EndNote
  • RefWorks
  • TxT
Links
Articles
Latest Articles
Most Read
Collections
Updates
Events
News
Multimedia
About
About Us
Contact
No. 86 Xueyuan South Road, Haidian District, Beijing
100081
010-62199257
qkjq@cast.org.cn

Copyright © 2025 China Association for Science and Technology. All rights reserved. For all open access content, the relevant licensing terms apply.
Sponsored by the Office of the Leading Group for Cybersecurity and Informatization of CAST, and supported by Science and Technology Review Publishing House