2.1 Equations of motion

In this paper, the earth is considered as a rotating spherical model. Then the three degree of freedom point mass equations of motion of an entry vehicle are as follows:

(1) \begin{align}\left\{ \begin{array}{l}\dfrac{{\textrm{d}r}}{{\textrm{d}t}} = V\sin \theta \\[8pt]\dfrac{{\textrm{d}\lambda }}{{\textrm{d}t}} = \dfrac{{V\cos \theta \sin {\psi _V}}}{{r\cos \phi }}\\[8pt]\dfrac{{\textrm{d}\phi }}{{\textrm{d}t}} = \dfrac{{V\cos \theta \cos {\psi _V}}}{r}\\[8pt]\dfrac{{\textrm{d}V}}{{\textrm{d}t}} = - D - \dfrac{{\sin \theta }}{{{r^2}}} + \omega _e^2r\cos \phi (\!\sin \theta \cos \phi - \cos \theta \cos {\psi _V}\sin \phi )\\[8pt]\dfrac{{\textrm{d}\theta }}{{\textrm{d}t}} = \dfrac{1}{V}\left[L\cos \sigma + \left({V^2} - \dfrac{1}{r}\right)\dfrac{{\cos \theta }}{r} + 2{\omega _e}V\cos \phi \sin {\psi _V} + \right.\\[8pt] \left. \qquad\;\; \omega _e^2r\cos \phi (\cos \theta \cos \phi + \sin \theta \sin \phi \cos {\psi _V})\right]\\[8pt]\dfrac{{d{\psi _V}}}{{dt}} = \dfrac{1}{V}\left[\dfrac{{L\sin \sigma }}{{\cos \theta }} + \dfrac{{{V^2}}}{r}\cos \theta \sin {\psi _V}\tan \phi\; + \right. \\[8pt] \left. \qquad\quad 2{\omega _e}V(\cos \phi \cos {\psi _V}\tan \theta - \sin \phi ) + \dfrac{{\omega _e^2r\sin {\psi _V}\sin \phi \cos \phi }}{{\cos \theta }}\right],\end{array} \right.{}\end{align}

where $r$ is the radius from the earth center to the vehicle, $\lambda $ and $\phi $ are the longitude and latitude, respectively, $V$ is the velocity of the vehicle, $\theta $ is the flight-path angle, ${\psi _V}$ is the heading angle measured from the north in clockwise direction, ${\omega _e}$ is the rotational angular velocity of the earth, $\sigma $ is the bank angle, $L$ and $D$ are the lift and drag accelerations, respectively. All variables are dimensionless. $L$ and $D$ can be calculated by:

(2) \begin{align}\left\{ \begin{array}{l}L = \dfrac{1}{{2m{g_0}}}\rho {\left( {{V_c}V} \right)^2}{C_L}\left( {\alpha ,{V_c}V} \right){S_{ref}}\\[10pt]D = \dfrac{1}{{2m{g_0}}}\rho {\left( {{V_c}V} \right)^2}{C_D}\left( {\alpha ,{V_c}V} \right){S_{ref}}\end{array}, \right.{}\end{align}

where $m$ is the mass of the vehicle, ${g_0}$ is the gravitational acceleration of sea level, ${V_c} = \sqrt {{g_0}R} $ is the dimensionless parameter of flight speed, ${C_L}$ and ${C_D}$ are the lift and drag coefficients, respectively, $\alpha $ is the angle-of-attack, ${S_{ref}}$ is the reference area of the vehicle, $\rho $ is the air density which is calculated by:

(3) \begin{align}\rho \left( H \right) = {\rho _0}{e^{ - \beta \cdot H}}{,}\end{align}

where $H$ is the altitude, ${\rho _0} = 1.225\,{\rm{kg/}}{{\rm{m}}^{\rm{3}}}$ is the air density at sea level and $\beta = 1/7110$ .

2.2 Trajectory constraints

Typical reentry process constraints include:

1. (1) Path constraints:

(4) \begin{align}\dot Q = {C_Q}{\rho ^{0.5}}{\left( {{V_c}V} \right)^{3.15}} \le {\dot Q_{\max }}{,}\end{align}

(5) \begin{align}q = \frac{1}{2}\rho {\left( {{V_c}V} \right)^2} \le {q_{\max }}{,}\end{align}

(6) \begin{align}{n_y} = \sqrt {{L^2} + {D^2}} \le {n_{\max}}{,}\end{align}

where $\dot Q$ , $q$ , ${n_y}$ are the stagnation point heating rate, dynamic pressure and total aerodynamic load, respectively. And ${\dot Q_{\max }}$ , ${q_{\max }}$ , ${n_{\max }}$ are the maximum limits of them, respectively. ${k_Q} = 0.00015$ is a heat transfer coefficient. These three constraints are rigid constraints and must be satisfied.

In the lifting entry trajectory, the flight-path angle $\theta $ is small and varies slowly. Setting $\cos \theta = 1$ and $\dot \theta = 0$ in the dynamics equations of flight-path angle and ignoring the earth rotation gives:

(7) \begin{align}L\cos \sigma + \left( {{V^2} - \frac{1}{r}} \right)\frac{1}{r} = 0.\end{align}

Since $\cos \sigma \le 1$ must be satisfied, we obtain:

(8) \begin{align}L + \left( {{V^2} - \frac{1}{r}} \right)\frac{1}{r} \ge 0,\end{align}

which is called the quasi-equilibrium glide condition (QEGC) [Reference Shen and Lu3]. It means that only when the H–V profile is located below this QEGC boundary, the flight-path angle can be controlled by adjusting the bank angle.

1. (2) Terminal constraints:

In this paper, the terminal constraints are given as follows.

(9) \begin{align}\begin{array}{l}\left| {{r_f} - r_f^ * } \right| \le {\varepsilon _r}, \left| {{\lambda _f} - \lambda _f^ * } \right| \le {\varepsilon _\lambda } \\[4pt] \left| {{\phi _f} - \phi _f^ * } \right| \le {\varepsilon _\phi }, \left| {{V_f} - V_f^ * } \right| \le {\varepsilon _V},\end{array}\end{align}

where $\left( {{r_f},{\lambda _f},{\phi _f},{V_f}} \right)$ are the terminal states of the vehicle, $\left( {r_f^ * ,\lambda _f^ * ,\phi _f^ * ,V_f^ * } \right)$ are the desired terminal states, $\left( {{\varepsilon _r},{\varepsilon _\lambda },{\varepsilon _\phi },{\varepsilon _V}} \right)$ are the corresponding permissible errors.

1. (3) Geographic constraints:

No-fly zones and waypoints are two kinds of geographic constraint that shall be satisfied in entry trajectory planning. The vehicle must avoid the no-fly zones and cross from the waypoints. The waypoint constraint can be expressed as a point $\left( {{\lambda _w},{\phi _w}} \right)$ that should be passed. The terminal constraints can also be regarded as waypoint constraints. The no-fly zone model is an infinitely high cylinder, whose center is $\left( {{\lambda _n},{\phi _n}} \right)$ with a radius of ${R_n}$ . The no-fly zone constraint can be expressed as that all the points $\left( {{\lambda _M},{\phi _M}} \right)$ in the entry trajectory satisfying:

(10) \begin{align}{S_n}\left( {{\lambda _M},{\phi _M},{\lambda _n},{\phi _n}} \right) \gt {R_n},\end{align}

where ${S_n}\left( {{\lambda _M},{\phi _M},{\lambda _n},{\phi _n}} \right) = {R_0}\arccos \left( {\sin {\phi _M}\sin {\phi _n} + \cos {\phi _M}\cos {\phi _n}\cos \left( {{\lambda _M} - {\lambda _n}} \right)} \right)$ . ${R_0}$ is the radius of the earth.

1. (4) Control constraints:

The control variables in entry trajectory planning are usually the angle-of-attack $\alpha $ and bank angle $\sigma $ . And they are constrained within a specific range due to the practical flight conditions:

(11) \begin{align}{\alpha _{\min }} \le \alpha \le {\alpha _{\max }}, {\sigma _{\min }} \le \sigma \le {\sigma _{\max }},\end{align}

where ${\alpha _{\min }}$ , ${\alpha _{\max }}$ are the minimum and maximum constrained values of $\alpha $ , respectively, and $ {\sigma _{\min }}$ , ${\sigma _{\max }}$ are the minimum and maximum constrained values of $\sigma $ , respectively.

3.1 Constraint conversion and reentry corridor

This paper employs altitude-versus-velocity (H-V) plane to design the entry corridor and determine the maximum and minimum boundaries allowed. Since the boundaries of the constraints such as overload and QEGC are obtained under a certain angle-of-attack, it is necessary to design a nominal angle-of-attack in advance according to the capability of the vehicle. In this paper, the nominal angle-of-attack is defined as a function of velocity:

(12) \begin{align}\alpha = \left\{ \begin{array}{l@{\quad}l}{\alpha _{\max }}& V \gt {V_1}\\[4pt]\dfrac{{{\alpha _{L/{D_{\max }}}} - {\alpha _{\max }}}}{{{V_2} - {V_1}}}\left( {V - {V_1}} \right) + {\alpha _{\max }} & {V_2} \le V \le {V_1}\\[4pt]{\alpha _{L/{D_{\max }}}} & V \lt {V_2},\end{array} \right.{}\end{align}

where ${\alpha _{\max }}$ the maximum angle-of-attack, ${\alpha _{L/{D_{\max }}}}$ is the angle-of-attack which has the largest lift-to-drag ratio, ${V_1}$ and ${V_2}$ are two specific velocities. The purpose of designing the angle-of-attack function in (12) is to decelerate rapidly in the early stage of reentry, and increase the altitude of the first heat rate peaking point as much as possible to meet the heating rate constraints. When the vehicle crosses the peaking point of the first heat rate, it then flies with the maximum lift-to-drag ratio, which increase the cross range in the gliding phase. The schematic diagram of the nominal angle-of-attack is shown in the Fig. 1.

Figure 1. Nominal angle-of-attack.

Then the path constraints (4)–(7) can be transformed into an entry corridor in H-V plane. The upper boundary of the entry corridor is defined by QEGC and the lower boundary is defined by heating rate, dynamic pressure and aerodynamic load constraints. The entry corridor based on H-V plane is described as follows:

(13) \begin{align}\left\{ \begin{array}{l}{H_{\min }}(V) = \max \left( {{H_{{{\dot Q}_{\max }}}}(V),{H_{{q_{\max }}}}(V),{H_{{n_{{y_{\max }}}}}}(V)} \right)\\[4pt]{H_{\max }}(V) = {H_{QEGC}}(V),\end{array} \right.{}\end{align}

where ${H_{\min }}(V)$ , ${H_{\max }}(V)$ are the upper boundary and lower boundary of the entry corridor, respectively, ${H_{{q_{\max }}}}(V)$ , ${H_{{n_y}_{\max }}}(V)$ , ${H_{{{\dot Q}_{\max }}}}(V)$ , ${H_{QEGC}}(V)$ represent the altitude constraints of the dynamic pressure, aerodynamic load, heating rate and QEGC under a certain velocity, respectively. From (2)–(6), (13), a mathematical expression of the basic reentry corridor boundary can be determined, as follows:

(14) \begin{align}\begin{array}{l}H \gt \dfrac{2}{\beta }\ln \left( {\dfrac{{{C_Q}{V^{3.15}}}}{{{{\dot Q}_{s\max }}}}} \right) = {H_{{{\dot Q}_{max}}}}(V)\\[10pt]H \gt \dfrac{1}{\beta }\ln \left( {\dfrac{{{\rho _0}{V^2}}}{{2{q_{\max }}}}} \right) = {H_{{q_{max}}}}(V)\\[10pt]H \gt \dfrac{1}{\beta }\ln \left( {\dfrac{{{\rho _0}{V^2}S\sqrt {C_D^2 + C_L^2} }}{{2{n_{y\max }}m{g_0}}}} \right) = {H_{{n_{y\max }}}}(V).\end{array}{}\end{align}

QEGC in (7) has no analytical solution and needs to be solved by numerical method. In this paper, a dichotomy method is used for any given velocity, and a sufficiently accurate ${H_{QEGC}}$ is obtained by iterating 10 times in the range of 25–80 km. QEGC is a soft constraint and the reentry corridor obtained by the above constraints is shown in Fig. 2.

Figure 2. Reentry corridor.

In order to convert the QEGC into the bank angle boundary, the constraint conditions (4)–(6) of the reentry corridor are brought into Equation (7), and the bank angle function is obtained with velocity as the independent variable.

(15) \begin{align}\begin{array}{l}{\left| \sigma \right|_{QEG\_\dot Q}}(V) = {\left| \sigma \right|_{QEG\_\dot Q}}\left( {{H_{{{\dot Q}_{max}}}}(V),V} \right)\\[4pt]{\left| \sigma \right|_{QEG\_q}}(V) = {\left| \sigma \right|_{QEGQ\_q}}\left( {{H_{{q_{max}}}}(V),V} \right)\\[4pt]{\left| \sigma \right|_{QEG\_n}}(V) = {\left| \sigma \right|_{QEGQ\_n}}\left( {{H_{{n_{y\max }}}}(V),V} \right){.}\end{array}\end{align}

Then the upper boundary of the bank angle in the gliding phase is:

(16) \begin{align}{\left| \sigma \right|_{\max }}(V) = \min \left( \begin{array}{l}{\left| \sigma \right|_{QEG\_\dot Q}}\left( {{H_{{{\dot Q}_{max}}}}(V),V} \right)\\[4pt]{\left| \sigma \right|_{QEGQ\_q}}\left( {{H_{{q_{max}}}}(V),V} \right)\\[4pt]{\left| \sigma \right|_{QEGQ\_n}}\left( {{H_{{n_{y\max }}}}(V),V} \right).\end{array} \right){}\end{align}

Let ${\sigma _{QEGC}} = 0$ be the lower boundary of the bank angle in the gliding phase, therefore, for any velocity $V$ in the gliding phase, as long as it satisfies

(17) \begin{align}{\sigma _{QEGC}} \le \sigma(V) \le {\sigma _{\max }}(V).\end{align}

Equation (17) can guarantee the vehicle’s flight trajectory will always remain within the reentry corridor. From (17), the bank angle is restricted during the reentry process, and the bank angle corridor is shown in Fig. 3.

Figure 3. Bank angle corridor.

3.2 Bank angle profile analysis

The core of the predictor corrector algorithm is to design the bank angle profile while latitude and longitude constraints are usually transformed into range constraints in designing longitudinal trajectory, Thus it is necessary to analyse the relationship between the bank angle profile with the reentry range, altitude, velocity. For reentry range $s$ there is:

(18) \begin{align}\dot s = \frac{{V\cos \theta }}{r}.\end{align}

From (18) and the fourth equation in (1), the derivative of range to velocity can be obtained:

(19) \begin{align}\frac{{ds}}{{dV}} = - \frac{{Vr\cos \theta }}{{\sin \theta + D{r^2}}}.\end{align}

Given that $\sin \theta \approx 0$ , $\cos \theta \approx 1$ , there is:

(20) \begin{align}\frac{{ds}}{{dV}} = - \frac{V}{{Dr}}.\end{align}

Analyse (20) in the height-velocity (H-V) plane, due to the dimensionless radius from the earth center to the vehicle $r \approx 1$ , for a certain speed, the higher the altitude, the smaller the resistance acceleration, the larger the range change rate, and the larger the total cross range. That is, for a given H-V curve with specific initial and terminal velocity, the higher the flying height, the larger the corresponding range.

It can be obtained by transforming the QEGC to (20):

(21) \begin{align}\frac{{ds}}{{dV}} = - \frac{{VL/D}}{L} = - \frac{{L/D\cos \sigma V}}{{1/r - V}}.\end{align}

Assuming a constant lift-to-drag ratio of the reentry vehicle, and the integral of (21) is available:

(22) \begin{align}s = \frac{{L/D}}{2}\cos \sigma \ln \frac{{1/r - V_f^2}}{{1/r - V_0^2}}.\end{align}

Equation (22) shows that when giving the initial and terminal velocity, the range will depend on the lift-to-drag ratio and the bank angle. The smaller the bank angle, the larger the cross range.

In addition, for QEGC (7), since $r \approx 1$ , there is:

(23) \begin{align}L\cos \sigma + {V^2} - 1 = 0.\end{align}

Since the lift acceleration $L$ decreases as the height increases, the corresponding equilibrium gliding bank angle decreases. That is, the smaller the value of the bank angle, the larger the height.

Based on the above analyses, the QEGC of the vehicle has the following rules: When the initial and terminal velocities are given, the smaller the bank angle is, the larger the flight altitude is, and the larger the range is. The bank angle, flight altitude and range are one-to-one correspondence. Larger ranges often correspond to larger terminal altitudes and smaller bank angles. However, when the given range constraint is large, and the terminal altitude constraint is small, or the given range constraint is small and the terminal altitude constraint is large, the un-match between the range constraint and the altitude constraint occurs.

Traditional predictor corrector guidance typically uses a constant or linear bank angle and gives terminal range constraints, while terminal altitude and velocity constraints are converted to energy constraints. When the terminal range constraint and the altitude constraint are corresponding or matched, a bank angle profile can be found to satisfy the range constraint while satisfying the altitude constraint. However, when the terminal range constraint doesn’t match the altitude constraint, i.e. the required bank angles of the two are inconsistent (for example, a larger range constraint requires a smaller bank angle, while a lower height constraint requires a larger bank angle). Since the predictor corrector algorithm uses the terminal trajectory error as the corrector amount, the designed bank angle profile can only satisfy the range constraint, and the terminal altitude can’t be accurately controlled. The algorithm for the mismatch of altitude and cross range is proposed below.

3.3 Longitudinal guidance

The initial descent phase has a high-flying altitude and a weak aerodynamic force. It is usually guided by a constant bank angle which is a kind of open-loop control, and transfers to the glide phase when certain conditions are met.

For the glide phase, based on the analyses of the previous section, the bank angle profile for the case where the terminal range constraint doesn’t match the altitude constraint in longitudinal plane is proposed in this section. In order to simplify the calculation and improve the efficiency of online calculation, the longitudinal motion equation with the dimensionless energy as the independent variable is established:

(24) \begin{align}\left\{ \begin{array}{l}\dfrac{{ds}}{{de}} = \dfrac{{\cos \theta }}{{rD}}\\[10pt]\dfrac{{dr}}{{de}} = \dfrac{{\sin \theta }}{D}\\[10pt]\dfrac{{d\theta }}{{de}} = \left[ {L\cos \sigma + \left( {{V^2} - \dfrac{1}{r}} \right)\dfrac{{\cos \theta }}{r}} \right]\dfrac{1}{{D{V^2}}},\end{array} \right.{}\end{align}

where $V = \sqrt {2\left( {1/r - e} \right)} $ . According to (24), the terminal range, altitude, and flight path angle can be predicted.

For the bank angle profile, in order to improve the ability to adjust the bank angle profile and meet the terminal altitude constraint, a segment linear profile is designed:

(25) \begin{align}\sigma = \left\{ \begin{array}{l}{\sigma _0} + \left( {e - {e_0}} \right)\dfrac{{{\sigma _0} - {\sigma _{mid}}}}{{{e_0} - {e_{mid}}}}, {e_0} \le e \lt {e_{mid}}\\[8pt]{\sigma _{mid}} + \left( {e - {e_{mid}}} \right)\dfrac{{{\sigma _f} - {\sigma _{mid}}}}{{{e_f} - {e_{mid}}}}, {e_{mid}} \le e \lt {e_f},\end{array} \right.{}\end{align}

where, ${e_0}$ , ${e_f}$ are the initial energy value and the terminal constraint energy value of the glide phase, respectively. ${e_{mid}} = \left( {{e_0} + {e_f}} \right)/2$ is the energy value at the midpoint, ${\sigma _0}$ , ${\sigma _f}$ , ${\sigma _{mid}}$ are the bank angle for the corresponding energy value. ${\sigma _0}$ and ${\sigma _f}$ are to be designed, ${\sigma _{mid}}$ is obtained by predictor corrector, and there is:

(26) \begin{align}e = \frac{1}{r} - \frac{{{V^2}}}{2}.\end{align}

Since the terminal altitude is closely related to the terminal bank angle, the above QEGC is used at the terminal point:

(27) \begin{align}L\left( {{V_f},{r_f}} \right)\cos {\sigma _f} + \left(V_f^2 - \frac{1}{{{r_f}}}\right)\frac{{\cos {\theta _f}}}{{{r_f}}} = 0.\end{align}

That is, the terminal bank angle ${\sigma _f}$ is selected as the QEGC bank angle under the terminal altitude and velocity constraint. The $\,{\sigma _f}\,$ is used as the initial value of the predictor corrector algorithm. In the guidance process, considering that the QEGC can’t be fully satisfied, in order to further improve the terminal accuracy, the dynamic equation of the flight path angle is obtained:

(28) \begin{align}{\dot \theta _f} = \frac{1}{{{V_f}}}\left[ {L\left( {{V_f},{r_f}} \right)\cos {\sigma _f} + \left(V_f^2 - \frac{1}{{{r_f}}}\right)\frac{{cos\;{\theta _f}}}{{{r_f}}}} \right],\end{align}

where ${\theta _f}$ is the terminal flight path angle predicted by (24) for the previous guidance period, And ${\dot \theta _f} = {{d{\theta _f}} \mathord{\left/ {\vphantom {{d{\theta _f}} {dt}}} \right.} {dt}}$ can also be predicted by Equation (24). Then, in each guidance period, it is continuously updated by (28), so that the terminal altitude constraint can be satisfied.

For ${\sigma _0}$ , first calculate a constant ${\sigma _s}$ that satisfies the range constraint. Comparing ${\sigma _s}$ with the bank angle ${\sigma _f}$ required for the terminal altitude constraints. Assume $\,{\sigma _s} \gt {\sigma _f}$ , in order to make the entire bank angle profile equal to that of the constant ${\sigma _s}$ profile which is for range constraints, take ${\sigma _s} \lt {\sigma _0} \lt {\sigma _{\max }}$ . While for $\,{\sigma _s} \lt {\sigma _f}$ , take $\,0 \lt {\sigma _0} \lt {\sigma _s}$ . So far, the bank angle profile has been designed, as shown in Fig. 4. The area of the fold line profile in Fig. 4 is equal to the constant value ${\sigma _s}$ profile, which means that the range constraint is satisfied; the terminal point is ${\sigma _f}$ , indicating that the terminal altitude constraint is satisfied.

Within each guidance period, given the initial $\,{\sigma _{mid}}$ , the terminal range is obtained by integrating the Equation (24) from the current state to the terminal state. Since ${\sigma _{mid}}$ determines the entire bank angle profile and the terminal range, and it can be updated through (29) (called secant method) to satisfy the range constraints:

(29) \begin{align}\sigma _{mid}^{\left( {i + 1} \right)} = \sigma _{mid}^{\left( i \right)} - \frac{{\sigma _{mid}^{\left( i \right)} - \sigma _{mid}^{\left( {i - 1} \right)}}}{{s_f^i - s_f^{\left( {i - 1} \right)}}}\left( {s_f^i - s_f^ * } \right)\!,\end{align}

where $s_f^i$ is the predicted terminal range, $s_f^ * $ is the desired range-to-go, and they can be calculated from the arc length of the sub-satellite point.

Figure 4. Bank angle profile.

3.4 Lateral guidance

The magnitude of $\sigma $ can determine the longitudinal states $\left( {H,V,\theta } \right)$ . Then the lateral states $\left( {\lambda ,\phi ,\psi } \right)$ are dominated by the sign of $\sigma $ . Therefore, we need to design a reversal logic of bank angle to control the lateral states so as to satisfy geographic constraints. In this paper, the reversal logic is designed on the base of dynamic heading angle corridor.

Define ${\psi _{LOS}}$ as the line of sight (LOS) of the vehicles’ current position to the target, and the equation is:

(30) \begin{align}{\psi _{_{LOS}}} = \arctan \frac{{\sin \left( {{\lambda _f} - \lambda } \right)}}{{\cos \phi \tan {\phi _f} - \sin \phi \cos \left( {{\lambda _f} - \lambda } \right)}}.\end{align}

As for no-fly zone constraints, first judge whether the trajectory of the vehicle to the reference point is affected by the no-fly zone. When the vehicle and the terminal point is on the same side of the no-fly zone, the trajectory is not affected by the no-fly zone. Taking the current longitude and latitude of the vehicle as the judging standards, it can be divided into four situations:

(31) \begin{align}\left\{ \begin{array}{l}{\lambda _M} \gt {\lambda _n} + {R_n}/{R_0}\\[4pt]{\lambda _M} \lt {\lambda _n} - {R_n}/{R_0}\\[4pt]{\phi _M} \gt {\lambda _n} + {R_n}/{R_0}\\[4pt]{\phi _M} \lt {\lambda _n} - {R_n}/{R_0},\end{array} \right.{}\end{align}

where $\,Z\left( {{\lambda _n},{\phi _n}} \right)$ is the central position of a no-fly zone with a radius of $\,{R_n}$ .

When the trajectory is affected by the no-fly zone, the relative positions between the vehicle and the reference point are as follows.

For a no-fly zone $\,Z\left( {{\lambda _n},{\phi _n}} \right)$ , there are two tangent lines of the no-fly zone that cross the point $M\left( {{\lambda _M},{\phi _M}} \right)$ , and two tangent angles ${\psi _{MZ\_\min }}$ , ${\psi _{MZ\_\max }}$ can be determined. ${\psi _{MZ\_\min }}$ is the smaller one and ${\psi _{MZ\_\max }}$ is the larger one.

If $\left| {{\psi _{MZ\_\min }} - {\psi _{LOS}}} \right| \ge \left| {{\psi _{MZ\_\max }} - {\psi _{LOS}}} \right|$ , the vehicle should fly over the lower side of the no-fly zone, which is shown as Fig. 5(a). In this case, the feasible heading angle corridor is $\left[ {{\psi _{\min }},{\psi _{\max }}} \right] = \left[ {{\psi _{MZ\_\max }},\pi } \right]$ . Similarly, if $\left| {{\psi _{MZ\_\min }} - {\psi _{LOS}}} \right| \lt \left| {{\psi _{MZ\_\max }} - {\psi _{LOS}}} \right|$ , the vehicle should fly over the upper side of the no-fly zone, which is shown as Fig. 5(b). And the feasible heading angle corridor is $\left[ {{\psi _{\min }},{\psi _{\max }}} \right] = \left[ {0,{\psi _{MZ\_\min }}} \right]$ .

Figure 5. Heading anglecorridor for no-fly zones.

Figure 4 depicts the relationship between the heading angle corridor defined by the no-fly zone and the heading angle corridor defined by the terminal point. Taking into account the influences of the no-fly zone and terminal point comprehensively. it is possible to give a heading angle corridor satisfying all no-fly zone constraints and terminal constraints:

If ${\lambda _M} \lt {\lambda _n} - {R_n}/{R_0}$ , then the no-fly zone should be considered.

(1) If $\left| {{\psi _{MZ\_\min }} - {\psi _{LOS}}} \right| \ge \left| {{\psi _{MZ\_\max }} - {\psi _{LOS}}} \right|$ , then:

(32) \begin{align}\left\{ \begin{array}{l}{\psi _{\min }} = {\psi _{MZ\_\max }}\\[4pt]{\psi _{\max }} = {\psi _{MZ\_\max }} + \delta \psi .\end{array} \right.{}\end{align}

If $\left| {{\psi _{MZ\_\min }} - {\psi _{LOS}}} \right| \lt \left| {{\psi _{MZ\_\max }} - {\psi _{LOS}}} \right|$ , then:

(33) \begin{align}\left\{ \begin{array}{l}{\psi _{\min }} = {\psi _{MZ\_\max }} - \delta \psi \\[4pt]{\psi _{\max }} = {\psi _{MZ\_\max }}.\end{array} \right.{}\end{align}

The sign of the bank angle is:

(34) \begin{align}{\rm{sign}}\left( {{\sigma _n}} \right) = \left\{ \begin{array}{l} - 1, \psi \gt {\psi _{\max }}\\[4pt]1, \psi \lt {\psi _{\min }}\\[4pt]{\rm{sign}}({\sigma _{n - 1}}), {\psi _{\min }} \lt \psi \lt {\psi _{\max }},\end{array} \right.{}\end{align}

where ${\sigma _n}$ is the bank angle of the current time and ${\sigma _{n - 1}}$ is the bank angle of the previous time.

If ${\lambda _M} \gt {\lambda _n} + {R_n}/{R_0}$ , then the terminal constraints should be considered.

Define the current heading angle error as $\psi - {\psi _{LOS}}$ . When the heading angle error exceeds the preset error corridor, the sign of the bank angle is changed; when the heading angle error does not exceed the error corridor, the sign of the tilt angle is kept unchanged.

(35) \begin{align}{\rm{sign}}\left( {{\sigma _n}} \right) = \left\{ \begin{array}{l} - 1, \psi - {\psi _{LOS}} \gt \delta \psi \\[4pt]1, \psi - {\psi _{LOS}} \lt - \delta \psi \\[4pt]{\rm{sign}}({\sigma _{n - 1}}), - \delta \psi \lt \psi - {\psi _{LOS}} \lt \delta \psi, \end{array} \right.{}\end{align}

where ${\sigma _n}$ is the bank angle of the current time and ${\sigma _{n - 1}}$ is the bank angle of the previous time. $\delta \psi $ is the preset heading angle error threshold.

The normal angle-of-attack has been given in advance. The lateral sub-planner gives the sign of the bank angle, while the longitudinal sub-planner offers the magnitude of the bank angle. Then the three-dimensional entry trajectory can be obtained by integrating the equations of motion (1).