Rapid trajectory planning of a reusable launch vehicle for airdrop with geographic constraints

Abstract

Online feasible trajectory generation for an airdrop unpowered reusable launch vehicle is addressed in this article. A rapid trajectory planning algorithm is proposed to satisfy not only the multiple path and terminal constraints but also the complex geographic constraints of waypoints and no-fly zones. Firstly, the lower and upper boundaries of the bank angle that implement all the path constraints are obtained based on the quasi-equilibrium glide condition. To determine the bank angle directly, a weighted interpolation of the boundaries is then developed, which provides an effective approach to simplify the planning process as a one-parameter search problem. Subsequently, three types of lateral planning algorithms are designed to determine the sign of the bank angle according to the requirements of waypoints passage, no-fly-zones avoidance, and terminal constraints in the airdrop process, and the convergence of these methods for passing over the waypoints and meeting the terminal conditions has been clarified and formally demonstrated. Considering the constraints in the actual airdrop flight missions, the planning trajectory is divided into several subphases to facilitate the application of corresponding algorithms. Finally, the performance of the proposed algorithm is assessed through three airdrop missions of reusable launch vehicle with different geographic constraints. Besides, the effectiveness of the algorithm is demonstrated by the Monte Carlo simulation results.

Introduction

In recent decades, owing to the advantages of quick support capability, ultra-long flight range, and high security in global strike and space transportation, the airdrop system of reusable launch vehicles (RLVs) has attracted widespread interest in both military and civilian applications. ¹ Compared with common vehicles, there are many distinctive characteristics for RLVs, especially the strong coupling and highly nonlinear dynamics, which make the modelling and control of RLVs a significant challenge. ² Meanwhile, the changing of the moment of inertia, center of gravity position, and aerodynamic moment caused by airdrop brings about uncertain dynamics. ³ Therefore, to significantly improve the flexibility and reliability of RLVs during the airdrop process, it is necessary and important to carry out accurate attitude and position control of flight in the presence of airdrop’s effect. ⁴

Current researches related to the airdrop mainly focused on the analysis of airdrop problem, ⁵ modeling of vehicle dynamics, ⁶ trajectory planning, ⁷ control law design, ⁸ simulation of airdrop process, ⁹ and so on. Gao et al. ⁵ studied the influence of perturbation on the airdrop system. Tao et al. employed an improved genetic algorithm to solve the homing trajectory planning problem for parafoil airdrop system. ⁷ Chen and Shi designed a backstepping and switch controller to guarantee the system stability and attitude tracking. ⁸ There are still many problems to be solved in the trajectory planning of the airdrop RLV except for satisfying various path and terminal constraints. For some special flight missions, the RLV is required to precisely overfly some predetermined positions on the ground which are referred as waypoints, such as telemetry stations that can validate the navigation of the vehicle, ¹⁰ enemy targets that need to be destroyed, or disaster-hit area where relief materials should be airdropped. ¹ In addition, some special regions on the RLV flight path should be kept away owing to geopolitical restrictions or security threats, such as large cities, fire blockade area, or military bases, which can be regarded as no-fly zones. ^11,12 Therefore, these waypoints and no-fly zones make the trajectory calculations for airdropping more complicated.

Numerous studies have been conducted to develop advanced trajectory planning algorithms for the airdrop RLV that can deal with multiple flight constraints, including path constraints, terminal constraints, and geographic constraints. The entry guidance for the Space Shuttle orbiter is a successful application of the drag acceleration profile generation that satisfies the specified entry terminal criteria while maintaining the flight within path constraints. ¹³ Later, some improvements have been developed. Two acceleration planning methods based on a reduced-order model are proposed by Mease et al., which consider both longitudinal and lateral motions of RLVs. ¹⁴ Saraf et al. developed the so-called evolved acceleration guidance logic for entry, which consists of two parts: a trajectory planner to generate reference drag acceleration and heading angle profiles and a tracking law to follow the two profiles simultaneously. ¹⁵ Moreover, to deal with the geographic constraints, a global discrete search strategy is presented by Yang and Zhao for the generation of four-dimensional trajectory of autonomous aerospace vehicles with obstacles, conflicts, and threats. ¹⁶ Then, Jorris and Cobb developed the autonomous trajectory planning techniques to satisfy waypoints and no-fly zone constraints simultaneously for the two-dimensional and three-dimensional (3-D) Common Aero Vehicle model. ^10,17 Besides, various optimization approaches, such as improved Gauss pseudospectral method ¹⁸ and multiphase convex programming, ¹⁹ have been applied to obtain the optimal reentry trajectory with the two geographic constraints for hypersonic vehicle. However, these trajectory planning/optimization algorithms for waypoint and no-fly-zone constraints are all conducted offline.

For the airdrop of RLV, it is crucial to rapidly generate a feasible trajectory to the airdrop positions online, while avoiding some potential threats, which can give full play to the quick support capability of RLV. Shen and Lu applied the quasi-equilibrium glide condition (QEGC) to an onboard trajectory design approach, which can generate a feasible 3-D flyable trajectory within several seconds, ²⁰ and filled a gap in the onboard suborbital entry trajectory planning. ²¹ Zhang et al. proposed a new approach to generate the longitudinal drag profile and lateral lift-to-drag profile simultaneously. ²² When taking the waypoint and no-fly zone constraints into account, the magnitude and sign of the bank angle are rapidly calculated in parallel based on a Newton–Raphson iteration scheme ²³ and also determined by two online sub-planners, respectively. ¹¹ A novel 3-D autonomous gliding guidance algorithm is presented by Guo et al., where two-bank-reversal logic and a dynamic heading-error corridor are designed to meet the geographic constraints. ²⁴ The 3-D analytical glide formulae are adopted to plan the reference trajectory of aerodynamic and inertial forces for no-fly-zone avoidance by Yu et al. ²⁵ These algorithms deal with the problem more comprehensively, and they all generate 3-D trajectories. However, the 3-D profiles generation needs a large amount of calculation, and the implementation process is quite complex. Besides, some intelligent methods are also applied to this area. Artificial potential field method is employed to avoid the no-fly zones by Zhang et al. ²⁶ Yong et al. applied the intelligence control method based on brain emotional learning to develop a predictor–corrector guidance algorithm, which can generate the reference trajectory onboard with waypoints. ²⁷ With the adjustment of angle of attack and nonreversal bank angle, the reference profile similar to Zhang is generated rapidly and can avoid no-fly zone effectively. ²⁸ These methods are flexible and easy to implement, but the maneuverability of the vehicle is not fully utilized, resulting in that only fewer constraints can be handled. Although many researches have been conducted on the rapid trajectory planning of RLVs for airdrop, there are few algorithms to analyze the problem of whether the computation time of the geographic constrained trajectory generation meets the requirement of online operation. For some of the above trajectory planning algorithms, whose computational complexity is relatively large, they are not fast enough to run onboard.

In this article, a new online trajectory planning approach is proposed to rapidly generate a feasible trajectory satisfying the path, terminal, waypoint, and no-fly zone constraints for the airdrop RLV. Unlike the existing method of converting the path constraints into the altitude or drag acceleration corridors, the QEGC is applied to describe these path constraints as the bank angle corridor in this article. Then, a weighted combination of the lower and upper boundaries of the bank angle provides an effective way to simplify the planning process as a one-parameter search problem, considering trajectory range and lateral maneuverability. Three types of lateral planning algorithms are developed according to the different requirements of waypoints passage, no-fly-zones avoidance, and terminal constraints satisfaction in the airdrop process. In addition, the convergence of the methods for passing over a waypoint and for meeting the terminal conditions is clarified and formally demonstrated. Finally, three nominal RLV airdrop fight missions with different geographic constraints are simulated, and to further evaluate the effectiveness of the planning algorithm, the Monte Carlo simulation is conducted.

The main contributions of this study are summarized as follows. First, the whole trajectory is divided into several subphases adaptively during the planning to reduce the integral range of each iteration and further expedite the calculation speed. The computation complexity and computation time of this algorithm satisfy the requirements of online operation, which is the major enhancement compared with the existing trajectory planning algorithms. Second, owing to the independence of the trajectory planning in each subphase, this proposed algorithm can deal with the airdrop flight with different geographic constraints or even without geographic constraints, which significantly improves the adaptability of the algorithm.

The rest of this article is organized as follows. The general trajectory planning problem considering multiple constraints is formulated in the second section. The third section describes the planning algorithms for waypoints passage, no-fly zone avoidance, and terminal constraints satisfaction. In the fourth section, simulation results of three nominal cases and some dispersion cases are presented to verify the effectiveness of the proposed planning algorithm. At last, some conclusions are made in the final section.

Problem formulation

In this section, the dynamics and kinematics characteristics of the RLV are described, and the trajectory planning problem with multiple path constraints, geographic constraints, and terminal constraints is formulated.

Dynamics

The three-degree of freedom (3-DOF) point-mass dynamics of the airdrop RLV over a spherical rotating Earth are described by the following equations of motion ^20,29

$\{\begin{array}{l}\frac{\text{d}r}{\text{d}t}=v\text{sin}\gamma \hfill \\ \frac{\text{d}\theta }{\text{d}t}=\frac{v\text{cos}\gamma \text{sin}\psi }{r\text{cos}\varphi }\hfill \\ \frac{\text{d}\varphi }{\text{d}t}=\frac{v\text{cos}\gamma \text{cos}\psi }{r}\hfill \\ \frac{\text{d}v}{\text{d}t}={{\omega }_e}^{2}r\text{cos}\varphi (\text{sin}\gamma \text{cos}\varphi -\text{cos}\gamma \text{sin}\varphi \text{cos}\psi )\hfill \\ -D-g\text{sin}\gamma \hfill \\ \frac{\text{d}\gamma }{\text{d}t}=\frac{L\text{cos}\sigma }{v}-\frac{g}{v}\text{cos}\gamma +\frac{v}{r}\text{cos}\gamma +2{\omega }_e\text{cos}\varphi \text{sin}\psi \hfill \\ +\frac{{\omega }_e^{2}r}{v}\text{cos}\varphi (\text{cos}\gamma \text{cos}\varphi +\text{sin}\gamma \text{cos}\psi \text{sin}\varphi )\hfill \\ \frac{\text{d}\psi }{\text{d}t}=\frac{L\text{sin}\sigma }{v\text{cos}\gamma }-2{\omega }_e(\text{tan}\gamma \text{cos}\psi \text{cos}\varphi -\text{sin}\varphi )\hfill \\ +\frac{v}{r}\text{cos}\gamma \text{sin}\psi \text{tan}\varphi +\frac{{\omega }_e^{2}r}{v\text{cos}\gamma }\text{sin}\psi \text{sin}\varphi \text{cos}\varphi \hfill \end{array}$

1

where the velocity heading angle ψ is measured from the North in a clockwise direction. L and D are described as

$\{\begin{array}{l}L=\frac{1}{2m}\rho v^2{S}_{\text{ref}}{C}_L\left(\alpha ,M_a\right),\\ D=\frac{1}{2m}\rho v^2{S}_{\text{ref}}{C}_D\left(\alpha ,M_a\right)\end{array}$

2

where C_L and C_D are determined by the angle of attack α and the Mach number M_a. The above notations are defined in Table 1. The geometric sketch of the airdrop flight is illustrated in Figure 1. The control variables are α and σ. To maintain the stability of RLV, it is necessary to limit the magnitude and rate of the two control variables.

Table 1.

Definition of common notations.

Notations	Definition and units
r	Radial distance from the Earth’s center to the vehicle, m
θ, ϕ	Longitude and latitude, respectively, deg
v	Earth-relative velocity, m/s
γ, ψ	Velocity flight-path angle and heading angle, respectively, deg
α, σ	Angle of attack and bank angle, respectively, deg
ωe	Self-rotation angular speed of Earth, rad/s
g	Gravitational acceleration, $g=\mu /r^2$ , ${\text{m/s}}^{\text{2}}$
μ	Gravitational constant, ${\text{m}}^{\text{3}}{\text{/s}}^{\text{2}}$
L, D	Aerodynamic lift and drag acceleration, respectively, ${\text{m/s}}^{\text{2}}$
CL, CD	Lift and drag coefficients, respectively
m	Vehicle mass, kg
$S_{ref}$	Reference area, ${\text{m/s}}^{\text{2}}$
re	Radius of Earth, m
h	Vehicle altitude from the sea level, $h=r-r_e$ , m
ρ	Density of atmosphere, $\rho ={\rho }_0\text{exp}(-h/h_s)$ , ${\text{kg/m}}^{\text{3}}$
ρ 0	Density of atmosphere at sea level, ${\text{kg/m}}^{\text{3}}$
hs	Scalar altitude coefficient
Ma	Mach number, $M_a=v/V_s$
Vs	Velocity of sound, m/s

Figure 1.

Geometric sketch of airdrop flight.

Path constraints

The typical inequality path constraints in the airdrop trajectory, consisting of the heating rate $\dot{Q}$ , dynamic pressure q, and aerodynamic load n, ¹¹ are given as

$\dot{Q}=c\sqrt{\rho }{v}^{3.15}\le {\dot{Q}}_{\text{max}}$

3

$q=\rho v^2/2\le q_{\text{max}}$

4

$n=\sqrt{L^2+D^2}\le n_{\text{max}}$

5

where c is the heating rate normalization constant. The parameters ${\dot{Q}}_{\text{max}}$ , $q_{\text{max}}$ , and $n_{\text{max}}$ represent the maximum allowable values of the heating rate, dynamic pressure, and aerodynamic load, respectively.

It is well-known that in a major portion of the flight trajectory, the flight path angle is small and varies relatively slowly. ²⁰ By setting $\gamma =0$ and $\dot{\gamma }=0$ and ignoring the Earth’s self-rotation in the fifth term of equation (1), a soft constraint, which need not to be strictly satisfied during the whole trajectory, is generated as

$L\text{cos}\sigma +\frac{v^2}{r}-g=0$

6

This is the so-called QEGC. It can guarantee that the airdrop trajectory doesn’t oscillate acutely for convenience of the attitude controller design and form the upper boundary of altitude corridor. Meanwhile, the above three inequality path constraints form the lower boundary of altitude corridor.

Geographic constraints

The geographic constraints generally include both waypoints and no-fly zone constraints for the airdrop process. Waypoints are defined as the predetermined positions on the ground for multiple enemy targets destruction or relief materials airdrop and are required to be precisely overflew by the airdrop RLV. ¹¹ That is, only the longitude and latitude for waypoint passage are constrained, and the other variables are uncertain. If the point $\left({\theta }_A,{\varphi }_A\right)$ on the flight trajectory is the closest point to the waypoint $\left({\theta }_{wp},{\varphi }_{wp}\right)$ , the constraint for the waypoint can be expressed as follows

${\left({\theta }_A-{\theta }_{wp}\right)}^2+{\left({\varphi }_A-{\varphi }_{wp}\right)}^2<{\epsilon }_{wp}^2$

7

where ${\epsilon }_{wp}$ is a preselected small positive value.

As mentioned before, the no-fly zones are specified as some special regions on the flight path that should not be violated owing to geopolitical restrictions or security threats. In this article, the no-fly zone is considered to be a cylinder zone with infinite altitude. ¹² If the distance from each position of the trajectory to the center of the no-fly zone is larger than the radius of the no-fly zone, the no-fly zone constraint is considered to be satisfied. ²⁴ Thus, the no-fly zone constraint can be given as follows

${\left(\theta -{\theta }_{nfz}\right)}^2+{\left(\varphi -{\varphi }_{nfz}\right)}^2>R_{nfz}^2$

8

where $\left(\theta ,\varphi \right)$ is the point on the flight trajectory, and $\left({\theta }_{nfz},{\varphi }_{nfz}\right)$ and $R_{nfz}$ are the predetermined center and radius of a no-fly zone, respectively. Note that the unit of $R_{nfz}$ is radians.

Terminal constraints

To ensure the successful flight in the terminal area energy management (TAEM) phase, the trajectory should satisfy some appropriate conditions at the TAEM interface. ^30,31 The typical terminal conditions on the altitude, velocity, and range to the target for the RLV airdrop flight are set as follows

$h_f=H_{\text{TAEM}},v_f=V_{\text{TAEM}},S_{\text{togo}\_f}=S_{\text{TAEM}}$

9

where h_f and v_f are the terminal altitude and velocity of RLV, respectively, and $S_{\text{togo}\_f}$ is the terminal range-to-go at the TAEM interface to the target auto landing interface (ALI). $H_{\text{TAEM}}$ , $V_{\text{TAEM}}$ , and $S_{\text{TAEM}}$ are all prescribed values of expectation.

In addition, it is desirable that the velocity vector at the TAEM interface should point closely toward the ALI. The terminal heading angle ψ_f constraint is given as

$|\Delta {\psi }_f|<{\epsilon }_{\psi }$

10

where $\Delta {\psi }_f={\psi }_f-{\psi }_{\text{TAEM}}$ . ${\psi }_{\text{TAEM}}$ denotes the desired heading angle from the final position to ALI, and ${\epsilon }_{\psi }$ is a specified tolerance.

Finally, too large a final bank angle $|{\sigma }_f|$ may lead to a long transient response for TAEM guidance and control. ³¹ Thus, the terminal bank angle constraint is defined as

$|{\sigma }_f|\le {\sigma }_{\text{TAEM}}$

11

where ${\sigma }_{\text{TAEM}}>0$ is a specified value.

The objective of the trajectory planning problem is to rapidly design a feasible trajectory from the initial position to the expected target, satisfying multiple flight constraints, including path constraints, terminal constraints, and geographic constraints for the airdrop RLV. It can be formulated as a feasible control problem, in which dynamics model is represented by equation (1), and multiple constraints are given in equations (3) and (11). Then, a novel online trajectory planning algorithm is developed in the next section to solve this problem.

Planning algorithm description

The total airdrop flight comprises the initial descent phase and the quasi-equilibrium glide (QEG) phase. ²⁰ During the initial descent phase, there is essentially no capability for carrying out large lateral maneuvers and the path constraints are not taken into account, due to the low atmospheric density. Starting at the initial interface, the RLV is required to descend and enter the airdrop flight corridor for transiting smoothly to a QEGC profile. This transition point represents the terminal of the initial descent phase. To generate the initial descent trajectory, the flexible planning algorithms have been studied ^11,20 and will not be discussed in this article. The QEG phase should ensure the satisfaction of the path constraints and terminal requirements, as well as executing lateral maneuvers to overfly the waypoints and avoid the no-fly zones. Therefore, we will focus on the trajectory planning of the glide phase below.

In this article, the nominal angle of attack profile is selected as piecewise linear functions of velocity and is described as follows

$\alpha \left(v\right)=\{\begin{array}{ll}{\alpha }_1,\hfill & v>V\hfill \\ \frac{{\alpha }_2-{\alpha }_1}{V_2-V_1}\left(v-V_1\right)+{\alpha }_1,\hfill & V_2\le v\le V_1\hfill \\ \frac{{\alpha }_f-{\alpha }_2}{V_{\text{TAEM}}-V_2}\left(v-V_2\right)+{\alpha }_2,\hfill & v<V_2\hfill \end{array}$

12

where V₁ and V₂ are the predefined velocity points of the profile, α₁ is the angle of attack preselected according to the heating rate constraint, while α₂ is predefined based on the range. The endpoint of the profile, α_f, is designed according to the QEGC with the designed final bank angle σ_f and the terminal conditions $V_{\text{TAEM}}$ and $H_{\text{TAEM}}$ .

Based on the specified angle of attack, the path constraints of equations (3) to (5) are converted into altitude constraints, with the minimum altitude boundary obtained as follows

$r_{\text{min}}\left(v\right)=\text{max}\left(r_{\dot{Q}}\left(v\right),r_q\left(v\right),r_n\left(v\right)\right)$

13

where, $r_{\dot{Q}}\left(v\right)$ , $r_q\left(v\right)$ , and $r_n\left(v\right)$ represent the altitude constraints corresponding to the maximum heating rate, dynamic pressure, and aerodynamic load, respectively. For a given value of v at a point on the trajectory, the QEGC defines the maximum allowable bank-angle as equation (14), based on the minimum altitude $r_{\text{min}}$

${\sigma }_{\text{QEGC}}\left(v\right)={{\displaystyle \text{cos}}}^{-1}\left(\frac{1}{L}\left(g-\frac{v^2}{r_{\text{min}}\left(v\right)}\right)\right)$

14

The QEGC provides a simple, yet effective method of limiting the bank-angle profile to implement all the path constraints. ${\sigma }_{\text{QEGC}}$ can be seen in Figure 2 as the main part of the upper boundary for $|\sigma |$ . To enforce the final bank angle constraint in equation (11), an appropriate threshold velocity V_t near the TAEM interface is chosen, and the upper boundary of bank angle is also parameterized as a linear function of velocity when the velocity is below V_t. From this, the velocity-dependent maximum allowable value of the bank angle, during the QEG phase, is expressed as

Figure 2.

The upper and lower boundaries of bank angle.

${\sigma }_{\text{up}}\left(v\right)\text{=}\{\begin{array}{l}{\sigma }_{\text{QEGC}}\left(v\right),V_t\le v\le V_{TP}\\ \left({\sigma }_{\text{QEGC}}\left(V_t\right)-{\sigma }_{up\_f}\right)\frac{\left(v-V_{\text{TAEM}}\right)}{V_t-V_{\text{TAEM}}}\\ +{\sigma }_{\text{up}\_f},V_{\text{TAEM}}\le v\le V_t\end{array}$

15

where V_TP is the velocity at the terminal of the initial descent phase. ${\sigma }_{\text{up}\_f}$ is the maximum final bank angle and should be set less than ${\sigma }_{\text{TAEM}}$ . This is expected to maintain the bank angle margin so that $|\sigma |\ge {\sigma }_{\text{down}}$ for most of the flight, which allows sufficient trajectory control to account for dispersions, where ${\sigma }_{\text{down}}<{\sigma }_{\text{TAEM}}$ is a preselected constant. Thus, the full constraints on the admissible magnitude of bank angle are obtained using equation (16) and shown in Figure 2

${\sigma }_{\text{down}}\le |\sigma \left(v\right)|\le {\sigma }_{\text{up}}\left(v\right)$

16

To observe the path constraints, a weighted combination of the lower and upper bounds of the bank angle is used to generate the bank angle profile. Along the airdrop flight trajectory, for any given v, the bank angle used for the trajectory integration is obtained as

$|\sigma \left(v\right)|\text{=}{\sigma }_{\text{down}}+k_{\sigma }\cdot \left({\sigma }_{\text{up}}\left(v\right)-{\sigma }_{\text{down}}\right)$

17

where $k_{\sigma }\in \left[0,1\right]$ denotes the interpolation coefficient. $k_{\sigma }=0$ indicates the minimum bank angle profile, and $k_{\sigma }=1$ the maximum. In addition to the path and terminal constraints, the lateral maneuverability requirements for waypoint passage and no-fly zone avoidance are considered in the design of $k_{\sigma }$ . To facilitate the design of α_f in equation (12), ${\sigma }_{\text{up}\_f}$ is set as ${\sigma }_{\text{up}\_f}\text{=}{\sigma }_{\text{down}}$ , such that the final bank angle σ_f will not be affected by $k_{\sigma }$ and can always have the same magnitude, $|{\sigma }_f|\text{=}{\sigma }_{\text{up}\_f}\text{=}{\sigma }_{\text{down}}$ .

The flowchart of the main blocks of the online trajectory planning algorithm is shown in Figure 3. After the planning of the initial descent phase, we need to determine the type of constraint that the RLV will face in the next subphase to adopt the corresponding planning algorithm. As the distance from the current position of the RLV to the predetermined waypoints, and to the tangent points of no-fly zones, can all be computed, the constraint closest to the current position is selected as the next constraint to be planned for. Only when all the geographic constraints have been successfully satisfied, the planning algorithm for terminal constraints will be employed.

Figure 3.

Flowchart of the online trajectory planning algorithm.

Passage strategy for waypoints

For each waypoint passage, the vehicle is allowed to cross the waypoint at different headings, hence the calculation of the expected heading at the waypoint presented in Liang et al. ³² is not required, which will save time. A single-bank-reversal logic is able to minimize the position error at the waypoint. ¹¹ The bank-reversal strategy, with an appropriate longitude θ_T as the bank-reversal point, is expressed as follows

$\text{sign}\left(\sigma \right)=\{\begin{array}{ll}{N}_{Init\_wp};\hfill & \theta <{\theta }_T\hfill \\ -N_{Init\_wp};\hfill & \theta \ge {\theta }_T\hfill \end{array}$

18

where ${\theta }_T\in \left[{\theta }_{Init},{\theta }_{wp}\right]$ . ${\theta }_{Init}$ is the initial longitude of this subphase. $N_{Init\_wp}$ is the initial sign of the bank angle and determined as follows

$N_{Init\_wp}=-\text{sign}\left({\psi }_{Init}-{\psi }_{ref}\right)$

19

where ${\psi }_{Init}$ is the initial heading, and ${\psi }_{ref}$ is the reference heading from the initial position toward the waypoint, which can be calculated based on the spherical–triangle theory according to the positional information

$\begin{array}{l}{\psi }_{ref}\text{=}{{\displaystyle \text{tan}}}^{-1}[\text{sin}\left({\theta }_{wp}-{\theta }_{Init}\right)/(\text{cos}\left({\varphi }_{Init}\right)\text{tan}\left({\varphi }_{wp}\right)\\ -\text{sin}\left({\varphi }_{Init}\right)\text{cos}\left({\theta }_{wp}-{\theta }_{Init}\right))]\end{array}$

20

With the magnitude and sign of the bank angle specified by equations (17) and (18), the dynamics equations (1) are numerically integrated from the current condition to ${\theta }_{wp}$ . The predicted final latitude ϕ_f is then compared to ${\varphi }_{wp}$ , and the mismatch is shown as follows

$f_{wp}={\varphi }_f-{\varphi }_{wp}$

21

Because θ_T has a large influence on the shape of the trajectory and determines the error f_wp, the planning process for waypoint passage is actually a root-finding problem of solving for θ_T to satisfy the condition $f_{wp}=0$ .

Theorem 1

The existence of θ_T to satisfy $f_{wp}=0$ is guaranteed, if the condition in equation (22) is satisfied

$f_{wp}({\theta }_{wp})\cdot N_{Init\_wp}<0$

22

where $f_{wp}\left({\theta }_{wp}\right)$ represents the predicted error f_wp based on the bank reversal point ${\theta }_T\text{=}{\theta }_{wp}$ .

Proof

When $N_{Init\_wp}=-1$ , according to the conditions in equations (18) and (19), we can obtain $f_{wp}({\theta }_{Init})<0$ , where $f_{wp}\left({\theta }_{Init}\right)$ represents the predicted error f_wp based on the bank reversal point ${\theta }_T\text{=}{\theta }_{Init}$ . It can be seen from the condition in equation (22) that $f_{wp}({\theta }_{wp})>0$ . Then, $f_{wp}({\theta }_{Init})\cdot f_{wp}({\theta }_{wp})<0$ . When $N_{Init\_wp}=1$ , according to the conditions in equations (18) and (19), we can get $f_{wp}({\theta }_{Init})>0$ . It can be seen from the condition in equation (22) that $f_{wp}({\theta }_{wp})<0$ . Then, $f_{wp}({\theta }_{Init})\cdot f_{wp}({\theta }_{wp})<0$ . Because f_wp is continuous in the domain ${\theta }_T\in \left[{\theta }_{Init},{\theta }_{wp}\right]$ and $f_{wp}({\theta }_{Init})\cdot f_{wp}({\theta }_{wp})<0$ is satisfied, the existence of θ_T to satisfy $f_{wp}=0$ can be guaranteed. This completes the proof of the theorem.

Two final longitudes, ${{\varphi }_f}^1$ and ${{\varphi }_f}^2$ , can be predicted according to the bank reversal points, ${{\theta }_T}^1$ and ${{\theta }_T}^2$ , respectively, where ${\theta }_{Init}\le {{\theta }_T}^1<{{\theta }_T}^2\le {\theta }_{wp}$ . When $N_{Init\_wp}=-1$ , as the schematic diagram shows in Figure 4, it can be easily understood that ${{\varphi }_f}^1<{{\varphi }_f}^2$ , which means $f_{wp}({{\theta }_T}^1)<f_{wp}({{\theta }_T}^2)$ . Thus, f_wp is monotonically increasing with θ_T(f_wp is monotonically decreasing with θ_T when $N_{Init\_wp}=1$ ) in the definition domain ${\theta }_T\in \left[{\theta }_{Init},{\theta }_{wp}\right]$ . The secant method below is effective in solving this root-finding problem, and the convergence of the algorithm for passing a waypoint is guaranteed

${{\theta }_T}^{\left(i+1\right)}={{\theta }_T}^{\left(i\right)}-\frac{{{\theta }_T}^{\left(i\right)}-{{\theta }_T}^{\left(i-1\right)}}{{f_{wp}}^{\left(i\right)}-{f_{wp}}^{\left(i-1\right)}}{f_{wp}}^{\left(i\right)}$

Figure 4.

Schematic diagram for waypoint passage.

where the superscripts $\left(i-1\right)$ , $\left(i\right)$ , and $\left(i+1\right)$ are the times of iteration.

Once this waypoint has been passed, the final states and control variables of this subphase can be regarded as the initial conditions of the next subphase. The same planning algorithm is then employed to pass the other waypoints, as per the flowchart in Figure 3.

Avoidance strategy for no-fly zones

As the bank reversal logic based on a reference point chosen near the no-fly zone, presented in Xie et al., ²³ cannot guarantee total avoidance effectively, the single-bank-reversal logic is unacceptable in this subphase. Furthermore, the reversal strategy should minimize the energy consumption caused by avoidance, so as to arrive at the terminal target successfully. ³³ Therefore, with the extension and simplification of the shuttle lateral guidance strategy, an avoidance strategy based on a dynamic heading-error corridor is used in this article. The most critical part of this strategy is the determination of the dynamic heading-error corridor, according to the relative position of the no-fly zone, reference point, and vehicle.

As seen in the schematic of a no-fly zone shown in Figure 5, we set the center and radius of the no-fly circular zone as $N\left({\theta }_{nfz},{\varphi }_{nfz}\right)$ and $R_{nfz}$ , respectively, $C\left(\theta ,\varphi \right)$ is the position of the RLV, and $T\left({\theta }_T,{\varphi }_T\right)$ is the reference point in the next subphase, which may be a waypoint or the target ALI. ψ_N and ψ_T are the headings from the current position of the RLV toward $N\left({\theta }_{nfz},{\varphi }_{nfz}\right)$ and $T\left({\theta }_T,{\varphi }_T\right)$ , respectively. ψ_A is the angle determined by the tangent from the RLV to the no-fly zone and influenced by the avoiding direction. All three angles are measured from the north in a clockwise direction and calculated in real time based on the current position of the RLV. The avoiding direction is indicated by an identifier Γ and defined in equation (23)

$\Gamma =-\text{sign}\left({\psi }_{Init}-{\psi }_N\right)$

23

Figure 5.

Schematic of no-fly zone.

where $\Gamma =1$ is the clockwise avoiding direction and $\Gamma =-1$ the counterclockwise avoiding direction.

For example, the avoiding direction in Figure 5 is counterclockwise according to the definition. It should be noted that the avoiding direction is determined independently at the beginning of each no-fly-zone avoidance subphase and will not change until this subphase is completed. On this basis, ψ_A is determined by the left tangent when $\Gamma =1$ and by the right tangent when $\Gamma =-1$ . The angle ψ_A can then be given by equation (24), uniformly

${\psi }_A={\psi }_N-\Gamma \cdot {\phi }_{\angle NCA}$

24

where ${\phi }_{\angle NCA}$ is the acute angle between $CA$ and $CN$ and calculated by equation (25)

${\phi }_{\angle NCA}={\text{sin}}^{-1}\left(\frac{R_{nfz}}{{\phi }_{CN}}\right)$

25

where ${\phi }_{CN}$ is the corresponding radian of the circle distance $CN$ .

$\begin{array}{l}{\phi }_{CN}={{\displaystyle \text{cos}}}^{-1}(\text{cos}\left({\varphi }_{nfz}\right)\cdot \text{cos}\left(\varphi \right)\cdot \text{cos}\left({\theta }_{nfz}-\theta \right)\\ +\text{sin}\left({\varphi }_{nfz}\right)\cdot \text{sin}\left(\varphi \right))\end{array}$

26

Therefore, the minimum and maximum boundaries of heading ψ can be obtained according to the avoiding direction. If $\Gamma \text{=}1$ , then

$\{\begin{array}{l}{\psi }_{\text{max}}=\text{min}\left({\psi }_A,{\psi }_T+\frac{1}{2}{\delta }_{\psi }\right)\hfill \\ {\psi }_{\text{min}}={\psi }_{\text{max}}-{\delta }_{\psi }\hfill \end{array}$

27

Else if $\Gamma \text{=}-1$ , then

$\{\begin{array}{l}{\psi }_{\text{min}}=\text{max}\left({\psi }_A,{\psi }_T-\frac{1}{2}{\delta }_{\psi }\right)\hfill \\ {\psi }_{\text{max}}={\psi }_{\text{min}}+{\delta }_{\psi }\hfill \end{array}$

28

where ${\delta }_{\psi }$ is a preselected positive value, which represents the magnitude of the heading-error corridor. It can be seen that the heading-error corridor dynamically changes with the variation of ψ_A and ψ_T in real time. Once the dynamic heading-error corridor of ψ is determined, the bank-reversal logic can be specified as follows

$\text{sign}\left(\sigma \right)=\{\begin{array}{ll}+1;\hfill & \psi <{\psi }_{\text{min}}\hfill \\ \text{sign}\left({\sigma }_{pre}\right);\hfill & {\psi }_{\text{min}}\le \psi \le {\psi }_{\text{max}}\hfill \\ -1;\hfill & \psi >{\psi }_{\text{max}}\hfill \end{array}$

29

where ${\sigma }_{pre}$ is the previous bank-angle command. The third term of equation (29) means that the vehicle maintains the sign of the previous bank angle when ψ is in the corridor.

Similar to the passage strategy for waypoints, the bank angle is generated by equations (17) and (29) and then numerically integrates the dynamics equations (1) from the current position until both of the following conditions are satisfied

$\Gamma \left({\psi }_A-{\psi }_T\right)>{\epsilon }_{TA}$

30

$\Gamma \left({\psi }_A-\psi \right)>{\epsilon }_{CA}$

31

where, ${\epsilon }_{TA}$ and ${\epsilon }_{CA}$ are preselected positive values. Equation (30) guarantees that the position of the RLV is away from the no-fly zone, while equation (31) ensures that the heading of the RLV is away from the no-fly zone. The final states and control variables of this subphase are then used as the initial conditions for the next subphase. The same planning algorithm is again employed to avoid the no-fly zones, as shown in the flowchart of Figure 3.

Plan for terminal constraints satisfaction

As per the flowchart in Figure 3, the planning algorithm for terminal constraint satisfaction will be employed when all the geographic constraints have been successfully satisfied. The TAEM interface point could be any position on the TAEM interface circle in which the terminal range-to-ALI is $S_{TAEM}$ , and the final desired longitude and latitude are not specified in advance. It is then unacceptable to apply the conventional bank-reversal logic based on the heading-error corridor that aims to fly the vehicle to a specified target position. In such a situation, we consider applying the single-bank-reversal logic by selecting an appropriate velocity V_T as the bank reversal point, instead of the longitude.

${\sigma }_{\text{Note}}=\{\begin{array}{ll}{N}_{Init\_ef};\hfill & v>V_T\hfill \\ -N_{Init\_ef};\hfill & v\le V_T\hfill \end{array}$

32

where $N_{Init\_ef}=\text{sign}\left({\psi }_{Init}-{\psi }_T\right)$ . $V_{Init}$ is the initial velocity of the final subphase. $N_{Init\_ef}$ is the initial sign of the bank angle and determined as follows

$N_{Init\_ef}=\text{sign}\left({\psi }_{Init}-{\psi }_{ref}\right)$

33

where ${\psi }_{ref}$ is the reference heading from the initial position to the target ALI.

The bank angle can be obtained by equation (17) and (32), and then the dynamics equations (1) are integrated from the initial state of the final subphase to $V_{\text{TAEM}}$ . The RLV can then definitely achieve the terminal velocity condition $v_f=V_{\text{TAEM}}$ . Because the magnitude of σ_f is a known constant and will not be affected by the weighted coefficient $k_{\sigma }$ , the design of the final angle of attack α_f can guarantee the terminal altitude condition $h_f=H_{\text{TAEM}}$ . Therefore, of the four terminal constraints, only $S_{\text{togo}\_f}$ and ψ_f will potentially not match the expected value at TAEM according to this planning algorithm. Based on the predicted final position $\left({\theta }_f,{\varphi }_f\right)$ , the desired heading ${\psi }_{\text{TAEM}}$ and the final range to the target ALI $S_{\text{togo}\_f}$ can be calculated. The corresponding terminal errors are then shown as follows:

$\{\begin{array}{l}{f}_S=S_{togo\_f}-S_{\text{TAEM}}\\ f_{\psi }={\psi }_f-{\psi }_{\text{TAEM}}\end{array}$

Because $k_{\sigma }$ and V_T together determine the terminal errors f_S and $f_{\psi }$ , the planning process for the terminal constraints is actually a one-parameter search problem of solving $\left(k_{\sigma },V_T\right)$ to satisfy the condition $f_S=0$ and $f_{\psi }=0$ . To facilitate the trajectory planning, and to simplify the computational complexity, $k_{\sigma }$ is mainly used for the correction of f_S as it has a large influence on the trajectory length, while V_T is mainly used for the correction of $f_{\psi }$ as it has a large influence on the shape of the trajectory.

Theorem 2

The existence of $k_{\sigma }$ to satisfy $f_S=0$ is guaranteed, if the conditions in equation (34) are satisfied.

$\{\begin{array}{l}{f}_S\left(0\right)<0\hfill \\ f_S\left(1\right)>0\hfill \end{array}$

34

where $f_S\left(0\right)$ and $f_S\left(1\right)$ represent the error f_S predicted based on the weighted coefficient $k_{\sigma }\text{=}0$ and $k_{\sigma }\text{=}1$ , respectively.

Proof

According to the conditions in equation (34), we can get $f_S\left(0\right)\cdot f_S\left(1\right)<0$ . Because f_S is continuous in the domain $k_{\sigma }\in \left[0,1\right]$ and $f_S\left(0\right)\cdot f_S\left(1\right)<0$ is satisfied, the existence of $k_{\sigma }$ to satisfy $f_S=0$ can be guaranteed. This completes the proof of the theorem.

Theorem 3

The existence of V_T to satisfy $f_{\psi }=0$ is guaranteed, if the condition in equation (35) is satisfied.

$f_{\psi }\left(V_{Init}\right)\cdot N_{Init\_ef}<0$

35

where $f_{\psi }\left(V_{Init}\right)$ represents the predicted error $f_{\psi }$ , based on the bank reversal point $V_T\text{=}{V}_{Init}$ .

Proof

When $N_{Init\_ef}=1$ , according to the conditions in equations (32) and (33), we can obtain $f_{\psi }\left(V_{\text{TAEM}}\right)>0$ , where $f_{\psi }\left(V_{\text{TAEM}}\right)$ represents the predicted error $f_{\psi }$ , based on the bank reversal point $V_T\text{=}{V}_{\text{TAEM}}$ . It can be seen from the condition in equation (35) that $f_{\psi }\left(V_{Init}\right)<0$ . Then, $f_{\psi }\left(V_{Init}\right)\cdot f_{\psi }\left(V_{\text{TAEM}}\right)<0$ . When $N_{Init\_wp}=-1$ , according to the conditions in equations (32) and (33), we can get $f_{\psi }\left(V_{\text{TAEM}}\right)<0$ . It can be seen from the condition in equation (35) that $f_{\psi }\left(V_{Init}\right)>0$ . Then, $f_{\psi }\left(V_{Init}\right)\cdot f_{\psi }\left(V_{\text{TAEM}}\right)<0$ . Because $f_{\psi }$ is continuous in the domain $V_T\in \left[V_{Init},V_{\text{TAEM}}\right]$ , and $f_{\psi }\left(V_{Init}\right)\cdot f_{\psi }\left(V_{\text{TAEM}}\right)<0$ is satisfied, the existence of V_T to satisfy $f_{\psi }=0$ can be guaranteed. This completes the proof of the theorem.

It can be easily understood from the planning algorithm that f_S is monotonically increasing with $k_{\sigma }$ . Under a certain $k_{\sigma }$ , the monotonicity of $f_{\psi }$ in the definition domain $V_T\in \left[V_{Init},V_{TAEM}\right]$ can be demonstrated in the same way as discussed in “Passage strategy for waypoints” section. A secant method is proven to be effective in solving this root-finding problem, and the convergence of the algorithm for meeting the terminal conditions is guaranteed.

$\{\begin{array}{l}{k_{\sigma }}^{\left(i+1\right)}={k_{\sigma }}^{\left(i\right)}-\frac{{k_{\sigma }}^{\left(i\right)}-{k_{\sigma }}^{\left(i-1\right)}}{{f_S}^{\left(i\right)}-{f_S}^{\left(i-1\right)}}{f_S}^{\left(i\right)}\\ {V_T}^{\left(i+1\right)}={V_T}^{\left(i\right)}-\frac{{V_T}^{\left(i\right)}-{V_T}^{\left(i-1\right)}}{{f_{\psi }}^{\left(i\right)}-{f_{\psi }}^{\left(i-1\right)}}{f_{\psi }}^{\left(i\right)}\end{array}$

Simulation results

According to the combination of two types of geographic constraints: waypoint and no-fly zone, there are three kinds of situations in the airdrop flight: The first case is only waypoints constraints, the second one is only no-fly zone constraints, and the third one is that two types of constraints are both encountered. These three situations basically include all the possible problems in the actual flight process. The only difference is the number of geographic constraints considered for each problem. Therefore, three corresponding mission scenarios for the airdrop flight of one kind RLV ^20,31 with different geographic constraints are developed for testing the proposed trajectory planning algorithm in this section.

The first mission is for flight with two waypoints and one no-fly zone constraints. The waypoints are located at (41.4°, 37.4°) and (59.4°, 47.4°). The no-fly zone is specified with the center at (56.4°, 36.4°) and a radius of 8°. In mission 2, five waypoints are predetermined, where the three waypoints are close to the target at longitudes and latitudes of (60.4°, 50.5°), (62.4°, 48.5°), and (64.4°, 46.5°), and the remaining two are relatively far at (47.4°, 41.5°) and (49.4°, 39.5°). The trajectory should pass one far waypoint and one close waypoint during the flight. Mission 3 is for flight with two no-fly zones. The no-fly zones are specified as follows: one with center (51.4°, 48.5°) and radius 8° and the other with center (67.4°, 41.5°) and radius 4.7°. Two cases are evaluated in mission 3 with the same initial conditions, except for the initial positions.

In the simulation, all the missions are constrained by the same peak heating rate of 800 ${\text{kW/m}}^{\text{2}}$ , peak dynamic pressure limit of 12 kPa, and peak aerodynamic load limit of 2.5 g. The initial altitude, velocity, flight-path angle, and heading angle are 120 km, 7515 m/s, −1.12°, and 49.32°, respectively. The initial longitude and latitude are generally set to (0°, 0°), only the setting of mission 3 is different. The same nominal angle of attack profile is used for all the three missions, where, $V_1=6000$ $\text{m/s}$ , $V_2=1500$ $\text{m/s}$ , ${\alpha }_1=40°$ , and ${\alpha }_2=22°$ . The relevant parameters of the bank angle profile are set as $V_t=1200$ $\text{m/s}$ , ${\sigma }_{\text{up}\_f}={\sigma }_{\text{down}}=10°$ . The altitude, velocity, and range to ALI at the nominal TAEM interface are 28 km, 900 m/s, and 100 km, respectively, with the errors less than 1 km, 5 m/s, and 5 km. The longitude and latitude of ALI are (69.4°, 51.3°). The maximum allowable TAEM heading error ${\epsilon }_{\psi }$ regarding the ALI is set as 5°, and the maximum value of the TAEM bank angle magnitude ${\sigma }_{TAEM}$ is 30°. In addition, the maximum magnitude of bank angle 80° and the maximum rate of bank angle 10°/s should be satisfied during the entire airdrop flight process. These constraints are implemented in the 3-DOF simulations.

The simulation environment are presented as follows:

Computer CPU: Intel Core i3-2120 3.30 GHz

Computer RAM: 4 GB DDR3 1333 MHz

Operating System: Windows 7 32-bit

Simulation Software: Microsoft Visual C++ 6.0

Nominal cases

The three nominal missions without any dispersions or uncertainties are simulated for the airdrop flight of RLV firstly, which is used to evaluate the performance of the proposed algorithm. The ground track for mission 1 in Figure 6 illustrates that the generated trajectory passes the waypoints and avoids the no-fly zone successfully. For the ground tracks in mission 2, shown in Figure 7, the four trajectories are obtained by passing different waypoint combinations with the same initial conditions. All the four trajectories satisfy the requirements of waypoint passage, with the position errors translated into distance on the ground, all less than 0.5 km. The ground tracks for mission 3 are shown in Figure 8. Trajectory 1 starts from the initial point where both the longitude and latitude are zero degrees, whereas that of trajectory 2 is (−29.2°, 42°). Two trajectories are generated with both no-fly zones avoided. The zoom-in view near the target ALI of the above three figures shows that the terminal range-to-go constraints of the three missions are well satisfied. The altitude and bank angle profiles for the three missions are shown in Figures 9 to 11. The altitude profiles illustrate that all path constraints, including waypoint and no-fly zone constraints, have been satisfied, and the trajectory phugoid oscillations can be further suppressed by broadening the boundaries of the bank angle if necessary. ³¹ It can be seen from the bank angle profiles that only a small number of bank reversals are conducted, and the magnitude of all the bank angle profiles decrease to approximately 10° toward the end. Table 2 lists the TAEM conditions for all three missions, which are also satisfied accurately. In the simulation, all the trajectories of the three missions are generated in 2–4 s (faster than other algorithms), which illustrates that the proposed trajectory planning algorithm can be used online for RLV airdrop flight guidance and control.

Figure 6.

Ground track for mission 1.

Figure 7.

Altitude and bank angle profiles for mission 1.

Figure 8.

Ground track for mission 2.

Figure 9.

Altitude and bank angle profiles for mission 2.

Figure 10.

Ground track for mission 3.

Figure 11.

Altitude and bank angle profiles for mission 3.

Table 2.

Terminal conditions achieved for the three missions.

Case		$h_f(\text{km})$	$v_f(\text{m/s})$	$S_{\text{togo}\_f}(\text{km})$	$\Delta {\psi }_f$ (deg)
Mission 1		28.05	899.76	99.73	3.28
Mission 2	Trajectory 1	28.15	899.79	100.92	−2.99
	Trajectory 2	28.08	899.81	99.54	−2.27
	Trajectory 3	28.12	899.90	100.82	4.04
	Trajectory 4	27.99	899.86	100.08	−3.18
Mission 3	Trajectory 1	28.03	899.57	99.46	−1.74
	Trajectory 2	28.16	899.61	100.53	0.38

Dispersion testing

To further demonstrate the effectiveness and robustness of the planning algorithm proposed in this article, a dispersion study is conducted through a Monte Carlo simulation with a variety of common random dispersions in the airdrop flight. The perturbations of initial conditions, atmospheric density, vehicle aerodynamic coefficients, and mass properties are taken into account in the additional 1000 dispersion cases, and all subject to Gaussian distribution. The mean and $3\sigma$ values of these dispersions are given in Table 3. The accuracy requirements of the terminal states are the same as nominal cases.

Table 3.

Statistics of dispersions for Monte Carlo simulations.

Dispersion parameter	Mean	$3\sigma$
Initial altitude (km)	0.0	1
Initial velocity (m/s)	0.0	50
Initial flight path angle (deg)	0.0	1
Initial heading angle (deg)	0.0	5
Atmospheric density	0.0	$10\%$
Aerodynamic coefficients	0.0	$5\%$
Mass	0.0	$5\%$

The dispersion simulation results are similar for all three missions, thus, only results of mission 1 are presented here. The ground tracks for mission 1 in the 1000 dispersion cases are presented in Figure 12, which shows that the trajectories have passed the waypoints and avoided the no-fly zone successfully. The altitude and bank angle profiles are plotted in Figure 13. The altitude profiles with the boundaries of the three path constraints illustrate the observance of the path constraints. The bank angle profiles show the satisfaction of the terminal bank angle constraint. Figure 14 illustrates a portion of the 100 km circle around the ALI and the tolerable error band of $\pm \text{5}$ km. The endpoints of the trajectories plotted for the 1000 dispersion cases are all lie within the $\pm \text{5}$ km band, which means that all the generated trajectories can reach the TAEM interface with acceptable errors of range to the target ALI. The terminal altitude and heading errors for the 1000 dispersion cases are shown in Figure 15, with the dispersion run numbers plotted on the abscissa. It can be seen from Figure 15 that no case exceeds the desirable $\pm \text{1}$ km altitude error or $\pm \text{5}°$ heading-error limits. The terminal conditions and path constraints obtained by Monte Carlo simulations are presented in Table 4. This table gives a more intuitively display of that all terminal conditions meet the requirements of precision and path constraints are satisfied accurately. The total simulation time is 3864.4 s and the average time of generating a trajectory is less than 4 s. These simulation results validate the robustness and rapidity of the algorithm.

Figure 12.

Ground track of the 1000 dispersion cases.

Figure 13.

Altitude and bank angle profiles of the 1000 dispersion cases.

Figure 14.

Terminal position for the 1000 dispersion cases.

Figure 15.

Terminal altitude and heading errors for the 1000 dispersion cases.

Table 4.

Terminal conditions and path constraints achieved for Monte Carlo simulations with dispersions.

Conditions	Maximum	Mean	Minimum	Variance	Reference
hf (km)	28.4598	28.0144	27.8601	0.0647	$\text{28}\pm \text{1}$
vf (m/s)	900.2526	899.7667	899.5468	0.1258	$\text{900}\pm \text{1}$
$S_{togo\_f}$ (km)	100.2998	99.7321	98.2212	0.6229	$\text{100}\pm \text{2}$
$\Delta {\psi }_f$ (deg)	4.8076	−1.1921	−3.9075	1.7938	$\text{0}\pm \text{5}$
$\text{max}(\dot{Q})({\text{kW/m}}^{\text{2}})$	783.0902	777.7928	771.8868	1.6941	$\text{< 800}$
$\text{max}(q)$ (kPa)	10.5118	10.147	9.2614	0.1529	<12
$\text{max}(n)$ (g)	2.1308	1.6595	1.6173	0.0584	$<2.5$

Conclusions

The rapid trajectory planning approach of the RLV developed in this study is intended for the application to airdrop flights with waypoint and no-fly-zone constraints. The performance of the planning algorithm is assessed through three different missions, in a high-fidelity simulation environment. The simulation on nominal cases provides an experimental verification of the rapid trajectory planning principles and an assessment of how well the algorithm work in different mission scenarios. In addition, Monte Carlo simulation results demonstrate the effectiveness of the algorithm on enforcing path and geographic constraints and achieving desired terminal conditions with high precision, even in the presence of some complex compound deviations. The robustness and the calculation speed of this algorithm are better than other algorithms, which makes it suitable for online operation. Future study will be conducted on the further reduction of the algorithm complexity and the improvement of its accuracy and adaptability.