# Boost-Phase Trajectory Planning with the Nonregular Reachable Area Constraints

Authors: Yin Diao; Rongjun Mu; Yingzi Guan; Naigang Cui
Journal: International Journal of Aerospace Engineering
 (2022)

---

## Abstract

During emergency return or mission change, the boost-glide vehicle needs to meet the nonregular reachable area constraints (NRACs). Therefore, this paper introduced NRACs into boost-phase trajectory planning to extend the maneuvering range and enhance the mission adaptability of boost-glide vehicle. Firstly, the library of the reachable area boundaries was constructed with the fast computation method under the deviations of boost terminal states, and the polynomials were used to fit the boundary parameters. Secondly, the nonlinear mapping relationship between the reachable area boundary parameters and the deviations of boost terminal states was obtained by the deep neural networks (DNNs). Then, the new boundary parameters were obtained by the constraint transformation rules under NRACs and the separation window constraints and transformed into the constraints of the boost-phase terminal states by DNNs. Finally, the hp-adaptive pseudospectral method (hpPM) was adopted to complete the trajectory planning considering the path and terminal constraints. The simulation results showed that the proposed trajectory planning method considering NRACs had high trajectory planning accuracy and good deviation adaptability and exhibited excellent performance in the adjustment of reachable area. This study provides theoretical support for integrate mission decisions and trajectory planning of boost-glide vehicles.


---

## Body

## 1. Introduction

Proposed by X. S. Qian in 1948, boost-glide vehicles, characterized with high lift-to-drag ratio, high altitude and velocity, and large flight envelope, are capable of long-range maneuvering, defense penetration, and target changing [1]. Diversified launch modes and multistage boosters have greatly expanded the glider’s application scope. In particular, the air-based launch mode allows the glider to quickly obtain the initial states in a specified direction, which brings great convenience for mission expansion. To reach the designated mission area precisely, suitable initial conditions of the glider are required for the boost-phase trajectory planning. Compared with the gliding phase, the boost phase is characterized by shorter flight time and weaker mobility, and it can significantly affect the reachable area by changing the initial conditions of the glider. Therefore, boost-phase trajectory planning has become one of the critical technologies affecting mission reachability.Under the preset separation conditions, boost-glide vehicles only depend on the gliding phase trajectory planning to meet the landing point constraints. In the case of emergency return of lift-body spacecraft with a booster during the ascending phase or in the case of the terminal mission change of the boost-glide vehicle during the boost stage, the spacecraft or the vehicle needs to meet the temporary constraints of terminal reachable area, and the adjustment ability of the gliding phase may not be enough to meet the constraints. Therefore, the authors consider using the remaining fuel in the booster to replan the boost-phase trajectory, make the glider reach more appropriate separation conditions, and make full use of the maneuverability of different flight phases to meet the terminal constraints. At the same time, the boost-phase trajectory is affected by the initial state of the vehicle, the combined action of aerodynamic force and engine thrust, and the environmental disturbances and needs to satisfy various constraints such as process state quantity, control angle, angular rate, and the quantity of terminal state quantity. Therefore, boost-phase trajectory planning is a complex practical problem to be addressed.With the rapid development of trajectory optimization theory [2–4], trajectory planning methods are mainly divided into methods based on optimal control principle, optimization theory, and artificial intelligence technology. Among them, the methods based on optimization theory have advantages in flexibility, but the online computing performance under complex constraints needs to be improved; the methods based on artificial intelligence technology have advantages in efficiency and adaptability [5]. The multiconstraint trajectory planning problem can be better solved with the rapid development of optimization theory and artificial intelligence technology [6]. In literature [7], a full-phase optimization model of the boost-glide vehicle was established, and the furthest range trajectory planning problem was solved via segment optimization. In literature [8–10], the gliding-phase trajectory planning problem satisfying the range and path constraints was investigated. In literature [11], the orthogonal test method was used to determine the optimal initial glide conditions and effective results were obtained. However, the connection between the reachable area and the boost terminal states was not established. In literature [12], the solution of the reachable area was transformed into a series of longest ranges to the virtual target points, and the optimal control theory and quasi-equilibrium glide condition were adopted to turn the longest ranges into univariate iterative problems, so that high computational accuracy and convergence speed were obtained. However, the actual coverage requirements were not considered in the boost-phase trajectory planning.DNNs/heuristic-based trajectory optimization algorithms have been widely studied. In reference [13, 14], a two-step strategy based on the trajectory optimization method and DNN attitude controller was proposed, which could better solve the problems in real-time attitude and trajectory planning and control for hypersonic vehicles. The desensitized/multiobjective trajectory optimization method generated the optimal reentry trajectory. In contrast, DNNs were used to obtain the optimal control attitude command matching the states, which could not realize the feedback from the reachable area constraints to the boost terminal states. Moreover, the thrust was not considered in the optimization model, so it could not be applied to the powered boost-phase trajectory optimization. In reference [15], a stochastic trajectory optimization framework was proposed to solve chance constraints, which could be used for the reentry trajectory planning of hypersonic vehicles under noise-perturbed dynamics and probabilistic constraints. In reference [16], a hybrid optimization algorithm composed of an initial value generator and an internal solver was proposed, which showed improved the convergence and stability. However, the algorithm proposed could not quickly calculate and adjust the reachable area.Although heuristic algorithms such as the multiobjective genetic algorithm [17] and gray wolf algorithm [18] have been widely studied and developed, the pseudospectral method is still one of the mainstream algorithms applied in the field of ascending trajectory planning [19]. On the basis of the equivalence between the Karush-Kuhn-Tucker (KKT) conditions of the parametric nonlinear program problem (NLP) and the discrete Hamiltonian boundary value problem (HBVP) [20], the Gaussian pseudospectral method with higher accuracy and convergence speed has been widely used in the boost-phase trajectory planning [21, 22]. With ever-increasing solution scale of trajectory planning problem under complex constraints in recent years, the hp-adaptive discrete grid updating method has been integrated into the Gaussian/Radau pseudospectral method for trajectory planning [23].To extend the maneuvering range of the boost-glide vehicle and meet NRACs, parametric descriptions of the reachable area, reachable area constraint transformation, and the boost-phase trajectory planning were studied in this paper. The main contributions of the work in this paper are drawn from the following three aspects.(1)
The NRACs were introduced into the boost-phase trajectory planning via constraint transformation to extend the maneuvering range and enhance mission adaptability(2)
A parametric method was designed to describe the reachable area boundary of the boost-glider vehicle, and DNNs were adopted to obtain the mapping relationship between the reachable area boundary and the booster’s separation point states(3)
The energy equivalence principles were proposed to deal with the separation window constraints, and the hp-adaptive pseudospectral method was adopted to plan the trajectory under path and terminal constraintsThe remainder of this paper is organized as follows. Firstly, the problem formulation and methodological framework are established in Section2. Secondly, the reachable area boundary library is established by the reachable area fast calculation method under initial deviations of states and environmental perturbations. And the boundaries of the reachable area are parameterized by the polynomial in Section 3. Next, constraints transformation rules, relational mapping of deep neural networks, and energy equivalence principles are proposed to transform the nonregular reachable area constraints into the terminal state constraints of the boost phase in Section 4. Then, the hpPM is adopted to plan the boost-phase trajectory under the path and terminal constraints in Section 5. Finally, the simulation analyses of the proposed method are carried out under the nonregular reachable area constraints in Section 6.


## 2. Problem Formulation and Methodological Framework

### 2.1. Problem Formulation

The boost-glide vehicle has two flight phases: the boost phase and gliding phase. The mathematical motion model of the boost phase can be established in the launch coordinate system as follows [24].
(1)mv̇=GBP00+GV−CxqSMCyqSMCzqSM+mgr′rx+R0xy+R0yz+R0z+mgωeωeωexωeyωez−mωex2−ωe2ωexωeyωexωezωexωeyωey2−ωe2ωezωeyωexωezωezωeyωez2−ωe2x+R0xy+R0yz+R0z−m0‐2ωez2ωey2ωez0‐2ωex−2ωey2ωex0ẋẏżẋẏż′=vṁ=−PIspg0,where v denotes the velocity vector; x, y, and z denote the three components of position, respectively; GB and GV denote the transformation matrix from the body-fixed coordinate system and velocity coordinate system to the launch coordinate system, respectively; P denotes the engine’s thrust; Cx, Cy, and Cz denote the drag, lift, and lateral force coefficients, respectively; q and SM denote the dynamic pressure and pneumatic reference area, respectively; r and ωe (which can be decomposed into ωex, ωey, and ωez) denote the geocentric vector diameter and earth rotation rate, respectively; gr′ and gωe denote the gravitational components along the direction of the geocentric vector and the direction of the earth’s rotation, respectively; g0 is the gravitational acceleration at sea level; Isp is specific impulse; ṁ is engine consumption per second.The dynamics model of the gliding phase is given as follows [12].
(2)ṙ=Vsinγ,θ̇=Vcosγsinψrcosφ,φ̇=Vcosγcosψr,V̇=−Dm+grsinγ+gωsinγsinφ+cosφcosψcosγ+rωE2cosφcosφsinγ−cosγsinφcosψ,γ̇=LcosσmV+Vrcosγ+grVcosγ+gωVsinφcosγ−cosφcosψsinγ+2ωEsinψcosφ+ωE2rcosφVcosφcosγ+sinφcosψsinγ,ψ̇=LsinσmVcosγ+Vrcosγsinψtanφ−gωVcosγcosφsinψ−2ωEtanγcosψcosφ−sinφ+ωE2rVcosγsinψsinφcosφ,where r=Re+h; Re denotes the radius of the earth, and h is the height; L and D denote the aerodynamic lift and drag, respectively; θ and φ denote the latitude and longitude, respectively; V is the velocity; γ and ψ denote the flight path angle and heading angle, respectively.The state equations of the boost terminal point need to be introduced to ensure the continuity of the motion equation.(3)V=vx2+vy2+vz2,r=x+R0x2+y+R0y2+z+R0z2,φ=arcsinr⋅ωerωe,θ=arctanrdyrdx,γ=arctanvyvx,ψ=arctanvENUxvENUy,where rdx and rdy denote the geocentric distance components in the X and Y directions in the geocentric coordinate system, respectively; vENUx and vENUy denote the velocity vector components in the X and Y directions in the local Cartesian coordinate coordinate system (ENU), respectively.GD denotes the coordinate conversion matrix from the launch coordinate system to the geocentric coordinate system, GE denotes the coordinate conversion matrix from the geocentric coordinate system to ENU. GD, GE, GB, and GV can be represented as follows [24].
(4)GD=Ly−90°+A0Lxφ0Lz−90°−θ0,GE=Lz180°Ly90°Ly−φLzθ,GB=LxγaLyψaLzφa,GV=LxνLyσvLzθv,where LiΔ,i=x,y,z is the direction cosine matrix of Δ rotation about i axis; φa, ψa, and γa denote the pitch angle, yaw angle, and roll angle, respectively; θv, σv, and ν denote the trajectory inclination angle, trajectory deflection angle, and velocity inclination angle, respectively; θ0, φ0, and A0 denote the initial latitude, longitude, and launch azimuth, respectively.


### 2.2. Methodological Framework

The boost-phase trajectory planning method considering NRACs proposed in this paper mainly consists of six parts. ① Polynomial was selected to parameterize the reachable area boundary obtained by the fast calculation method. ② The constraint transformation rules were proposed to calculate the new reachable area boundary parameters under NRACs. ③ DNNs were trained for the relationship mapping between the boost terminal states and the boundary parameters. ④ The initial values of the boost terminal states were obtained by the trained DNNs and the new reachable area boundary parameters. ⑤ Energy equivalence principles were proposed for the optimal boost terminal states under the separation window constraints. ⑥ The hpPM was used to plan the boost-phase trajectory considering path and terminal constraints. Therefore, the methodological framework is shown in Figure1.Figure 1
Methodological framework.


## 2.1. Problem Formulation

The boost-glide vehicle has two flight phases: the boost phase and gliding phase. The mathematical motion model of the boost phase can be established in the launch coordinate system as follows [24].
(1)mv̇=GBP00+GV−CxqSMCyqSMCzqSM+mgr′rx+R0xy+R0yz+R0z+mgωeωeωexωeyωez−mωex2−ωe2ωexωeyωexωezωexωeyωey2−ωe2ωezωeyωexωezωezωeyωez2−ωe2x+R0xy+R0yz+R0z−m0‐2ωez2ωey2ωez0‐2ωex−2ωey2ωex0ẋẏżẋẏż′=vṁ=−PIspg0,where v denotes the velocity vector; x, y, and z denote the three components of position, respectively; GB and GV denote the transformation matrix from the body-fixed coordinate system and velocity coordinate system to the launch coordinate system, respectively; P denotes the engine’s thrust; Cx, Cy, and Cz denote the drag, lift, and lateral force coefficients, respectively; q and SM denote the dynamic pressure and pneumatic reference area, respectively; r and ωe (which can be decomposed into ωex, ωey, and ωez) denote the geocentric vector diameter and earth rotation rate, respectively; gr′ and gωe denote the gravitational components along the direction of the geocentric vector and the direction of the earth’s rotation, respectively; g0 is the gravitational acceleration at sea level; Isp is specific impulse; ṁ is engine consumption per second.The dynamics model of the gliding phase is given as follows [12].
(2)ṙ=Vsinγ,θ̇=Vcosγsinψrcosφ,φ̇=Vcosγcosψr,V̇=−Dm+grsinγ+gωsinγsinφ+cosφcosψcosγ+rωE2cosφcosφsinγ−cosγsinφcosψ,γ̇=LcosσmV+Vrcosγ+grVcosγ+gωVsinφcosγ−cosφcosψsinγ+2ωEsinψcosφ+ωE2rcosφVcosφcosγ+sinφcosψsinγ,ψ̇=LsinσmVcosγ+Vrcosγsinψtanφ−gωVcosγcosφsinψ−2ωEtanγcosψcosφ−sinφ+ωE2rVcosγsinψsinφcosφ,where r=Re+h; Re denotes the radius of the earth, and h is the height; L and D denote the aerodynamic lift and drag, respectively; θ and φ denote the latitude and longitude, respectively; V is the velocity; γ and ψ denote the flight path angle and heading angle, respectively.The state equations of the boost terminal point need to be introduced to ensure the continuity of the motion equation.(3)V=vx2+vy2+vz2,r=x+R0x2+y+R0y2+z+R0z2,φ=arcsinr⋅ωerωe,θ=arctanrdyrdx,γ=arctanvyvx,ψ=arctanvENUxvENUy,where rdx and rdy denote the geocentric distance components in the X and Y directions in the geocentric coordinate system, respectively; vENUx and vENUy denote the velocity vector components in the X and Y directions in the local Cartesian coordinate coordinate system (ENU), respectively.GD denotes the coordinate conversion matrix from the launch coordinate system to the geocentric coordinate system, GE denotes the coordinate conversion matrix from the geocentric coordinate system to ENU. GD, GE, GB, and GV can be represented as follows [24].
(4)GD=Ly−90°+A0Lxφ0Lz−90°−θ0,GE=Lz180°Ly90°Ly−φLzθ,GB=LxγaLyψaLzφa,GV=LxνLyσvLzθv,where LiΔ,i=x,y,z is the direction cosine matrix of Δ rotation about i axis; φa, ψa, and γa denote the pitch angle, yaw angle, and roll angle, respectively; θv, σv, and ν denote the trajectory inclination angle, trajectory deflection angle, and velocity inclination angle, respectively; θ0, φ0, and A0 denote the initial latitude, longitude, and launch azimuth, respectively.


## 2.2. Methodological Framework

The boost-phase trajectory planning method considering NRACs proposed in this paper mainly consists of six parts. ① Polynomial was selected to parameterize the reachable area boundary obtained by the fast calculation method. ② The constraint transformation rules were proposed to calculate the new reachable area boundary parameters under NRACs. ③ DNNs were trained for the relationship mapping between the boost terminal states and the boundary parameters. ④ The initial values of the boost terminal states were obtained by the trained DNNs and the new reachable area boundary parameters. ⑤ Energy equivalence principles were proposed for the optimal boost terminal states under the separation window constraints. ⑥ The hpPM was used to plan the boost-phase trajectory considering path and terminal constraints. Therefore, the methodological framework is shown in Figure1.Figure 1
Methodological framework.


## 3. Parametric Description of the Reachable Area

### 3.1. Fast Calculation Method of the Reachable Area

The problem of the reachable area calculation can be decomposed into the issues on the longest ranges from the separation point to the virtual target points, as shown in Figure2. VTii=1,2⋯n denotes the virtual target points; v0 is the initial velocity direction of the vehicle; Fii=1,2⋯n represents the longest ground projection points under the energy management conditions, i.e., the farthest arrival points.Figure 2
Schematic diagram of the vehicle’s reachable area.When solving the issues on the longest range in the reachable area solution, the earth rotation, oblateness, and gravitational perturbation can be ignored, and the dynamic equations can be simplified based on the quasi-equilibrium glide condition.(5)θ̇=Vsinψcosϕ,ϕ̇=Vcosψ,V̇=−1CL/CD1−V2cosσ,ψ̇=1V1−V2tanσ+V2sinψtanϕ,where CL and CD denote the aerodynamic lift coefficient and drag coefficient in the gliding phase, respectively.The subproblem can be transformed into the following optimal control form:(6)minJs=min−cosSf,ẋs=fxs,u,t,H=λθθ̇+λϕϕ̇+λVV̇+λψψ̇,λ̇=−∂H∂xs∂H∂u=0,where Sf denotes the range; Js denotes the performance index; xs=θ,ϕ,V,ψ′ denotes the states; f⋅ denotes the right end term of the system differential equations; λθ, λϕ, λV, and λψ represent the corresponding covariate variables, respectively; H is Hamiltonian. The pitch angle σ is the only control variable when the attack profile α=αV is given in segments.Since the Hamiltonian function does not consider time, the analytical form of the optimal control solution can be obtained by using the control equation and the costate Equation (12):
(7)η1=sinχη2=−η1tanΦTη3=cosχ,tanσ=1−1/V2η1sinϕ+η2cosθ+η3sinθcosϕη1cosϕsinψ+η2sinθcosψ−cosθsinψsinϕ−η3sinθsinϕsinψ+cosθcosψsinχsinϕf−sinχtanΦTcosθf−ΘTcosϕf+cosχsinθf−ΘTcosϕf=0,where ΘT and ΦT denote the latitude and longitude of the virtual target point, respectively; θf and ϕf denote the terminal latitude and longitude of the gliding phase, respectively. The optimal control solution satisfying the requirements of the longest range can be found by iterating the univariate χ. The reachable area boundary can be obtained by connecting the farthest arrival point obtained by solving the issues on the longest range.The nominal initial states and deviations are designed as shown in Tables1 and 2.Table 1
Nominal states.
ItemValue (unit)Height70 kmVelocity5000 m/sFlight path angle-1°Heading angle90°Table 2
List of initial state deviations.
ItemMaximum value (unit)Height deviation±20 kmVelocity deviation±200 m/sFlight path angle deviation±5°The reachable area shown in Figure3 can be obtained by the fast calculation method by superimposing the normal distribution deviation shown in Table 2 on the nominal initial states. At the same time, the coordinate data of the reachable area is normalized to facilitate the fitting calculation of large data samples.Figure 3
The reachable area under nominal conditions.


### 3.2. Parametric Description

Considering that the boundary of the reachable area is a continuous closed curve connected by a series of boundary points, it can be described by curve parameters. Suppose that the boundary of the reachable area can be defined in the form of the following curve:(8)y=a+b1x+b2x2+⋯+csinwx+de−ξx+ηe−x−δ1/δ22.The least-square method is used to select the best fitting curves in different forms. The fitting results of different forms of curves are shown in Table3.Table 3
List of fitting results of different forms of curves.
OrderCurve formR-squareRMSE1y=a+b1x+b2x20.935600.042512y=a+b1x+b2x2+b3x30.935600.044403y=a+b2x20.935600.040844y=a+b2x2+csinwx0.935600.046565y=a+b2x2+de−ξx0.935600.044406y=a+b2x2+ηe−x−δ1/δ220.932300.04552The fitting results show that too much increase in the order and complexity of the curve will increase the fitting deviation. Therefore, the curve form corresponding to No. 3 is finally selected as the parametric description form of the reachable area boundary. The results after parameter normalization are shown in Figure4.Figure 4
Fitting results of the reachable area under nominal conditions. (Note:x and y denote the normalized longitudinal and cross ranges, respectively).


## 3.1. Fast Calculation Method of the Reachable Area

The problem of the reachable area calculation can be decomposed into the issues on the longest ranges from the separation point to the virtual target points, as shown in Figure2. VTii=1,2⋯n denotes the virtual target points; v0 is the initial velocity direction of the vehicle; Fii=1,2⋯n represents the longest ground projection points under the energy management conditions, i.e., the farthest arrival points.Figure 2
Schematic diagram of the vehicle’s reachable area.When solving the issues on the longest range in the reachable area solution, the earth rotation, oblateness, and gravitational perturbation can be ignored, and the dynamic equations can be simplified based on the quasi-equilibrium glide condition.(5)θ̇=Vsinψcosϕ,ϕ̇=Vcosψ,V̇=−1CL/CD1−V2cosσ,ψ̇=1V1−V2tanσ+V2sinψtanϕ,where CL and CD denote the aerodynamic lift coefficient and drag coefficient in the gliding phase, respectively.The subproblem can be transformed into the following optimal control form:(6)minJs=min−cosSf,ẋs=fxs,u,t,H=λθθ̇+λϕϕ̇+λVV̇+λψψ̇,λ̇=−∂H∂xs∂H∂u=0,where Sf denotes the range; Js denotes the performance index; xs=θ,ϕ,V,ψ′ denotes the states; f⋅ denotes the right end term of the system differential equations; λθ, λϕ, λV, and λψ represent the corresponding covariate variables, respectively; H is Hamiltonian. The pitch angle σ is the only control variable when the attack profile α=αV is given in segments.Since the Hamiltonian function does not consider time, the analytical form of the optimal control solution can be obtained by using the control equation and the costate Equation (12):
(7)η1=sinχη2=−η1tanΦTη3=cosχ,tanσ=1−1/V2η1sinϕ+η2cosθ+η3sinθcosϕη1cosϕsinψ+η2sinθcosψ−cosθsinψsinϕ−η3sinθsinϕsinψ+cosθcosψsinχsinϕf−sinχtanΦTcosθf−ΘTcosϕf+cosχsinθf−ΘTcosϕf=0,where ΘT and ΦT denote the latitude and longitude of the virtual target point, respectively; θf and ϕf denote the terminal latitude and longitude of the gliding phase, respectively. The optimal control solution satisfying the requirements of the longest range can be found by iterating the univariate χ. The reachable area boundary can be obtained by connecting the farthest arrival point obtained by solving the issues on the longest range.The nominal initial states and deviations are designed as shown in Tables1 and 2.Table 1
Nominal states.
ItemValue (unit)Height70 kmVelocity5000 m/sFlight path angle-1°Heading angle90°Table 2
List of initial state deviations.
ItemMaximum value (unit)Height deviation±20 kmVelocity deviation±200 m/sFlight path angle deviation±5°The reachable area shown in Figure3 can be obtained by the fast calculation method by superimposing the normal distribution deviation shown in Table 2 on the nominal initial states. At the same time, the coordinate data of the reachable area is normalized to facilitate the fitting calculation of large data samples.Figure 3
The reachable area under nominal conditions.


## 3.2. Parametric Description

Considering that the boundary of the reachable area is a continuous closed curve connected by a series of boundary points, it can be described by curve parameters. Suppose that the boundary of the reachable area can be defined in the form of the following curve:(8)y=a+b1x+b2x2+⋯+csinwx+de−ξx+ηe−x−δ1/δ22.The least-square method is used to select the best fitting curves in different forms. The fitting results of different forms of curves are shown in Table3.Table 3
List of fitting results of different forms of curves.
OrderCurve formR-squareRMSE1y=a+b1x+b2x20.935600.042512y=a+b1x+b2x2+b3x30.935600.044403y=a+b2x20.935600.040844y=a+b2x2+csinwx0.935600.046565y=a+b2x2+de−ξx0.935600.044406y=a+b2x2+ηe−x−δ1/δ220.932300.04552The fitting results show that too much increase in the order and complexity of the curve will increase the fitting deviation. Therefore, the curve form corresponding to No. 3 is finally selected as the parametric description form of the reachable area boundary. The results after parameter normalization are shown in Figure4.Figure 4
Fitting results of the reachable area under nominal conditions. (Note:x and y denote the normalized longitudinal and cross ranges, respectively).


## 4. Reachable Area Constraint Transformation

### 4.1. Constraint Transformation Rules

This section describes the method to transform NRACs to generate new coverage envelopes about the reachable area. As shown in Figure5, the projection of the boost terminal point on the ground is the origin O, and the coordinate system is established in the horizontal plane, with the axis Ox in the same direction as the projection of the initial velocity direction v0 on the ground. The original terminal states of the boost phase are h0,v0,γ0; x=a0y2+b0−ymax≤y≤ymax is the reachable area boundary curve generated by h0,v0,γ0; the nonregular reachable area is Q; the distance from the farthest point in Q to the coordinate origin is R; the angle between lL and lR formed by the left and right boundary points of Q and the coordinate origin O is Δξ; the centerline of Δξ is the new axis Ox′; the angle between Ox′ and Ox is Δψt; LT:x=ay2+b denotes the new reachable area boundary curve, so the key of the issue is to solve the parameters of LT.Figure 5
Description of the coverage requirements of the nonregular reachable area.Whether the original terminal states of the boost phase are appropriate can be judged by the relative position relationship between the boundary points ofQ and the boundary of the original reachable area.
(9)Judgment1:a0RsinΔψt2+b0>RcosΔψt,Judgment2:a0l0sinΔψt−12Δξ+Δξ0n2+b0>l0cosΔψt−12Δξ+Δξ0n,Judgment3:a0l1sinΔψt−12Δξ+Δξ1n2+b0>l0cosΔψt−12Δξ+Δξ1n⋯Judgmenti:a0li−2sinΔψt−12Δξ+Δξi−2n2+b0>li−2cosΔψt−12Δξ+Δξi−2n,where n is the number of divisions when judging the boundary for Q; li is the distance from the coordinate origin to the i‐th boundary judgment point. The n+2 judgment conditions must be satisfied simultaneously to determine whether the original terminal states of the boost phase are appropriate. The number of divisions can be obtained according to the description accuracy of the reachable area and the shape of Q.When the original terminal states of the boost phase do not meet NRACs, the change of initial heading angle can bring great convenience to the adjustment of the reachable area. However, the reachable area adjusted by the initial flight direction may still not meet the coverage requirements. Therefore, the judgments need to be reconducted in the newOx′y′-coordinate system.
(10)Judgment1:b0>R,Judgment2:a0l0sin−12Δξ+Δξ0n2+b0>l0cos−12Δξ+Δξ0n,Judgment3:a0l1sin−12Δξ+Δξ1n2+b0>l0cos−12Δξ+Δξ1n,⋯Judgmenti:a0li−2sin−12Δξ+Δξi−2n2+b0>li−2cos−12Δξ+Δξi−2n.When then+2 judgment conditions cannot be satisfied simultaneously, the new boundary parameters can be adjusted and determined according to the range of values of the boundary parameters.
(11)amin≤a≤amax,bmin≤b≤bmax,ymaxmin≤ymax≤ymaxmax.As a supporting parameter, the parameterymax varies inversely and equivalently with the adjustment of a. The adjustment strategies for the three parameters are shown in Figure 6.Figure 6
Adjustment strategies for the reachable area boundary parameters.


### 4.2. Relational Mapping of Deep Neural Networks

The relationship between the reachable area boundary parameters and the boost terminal states needs to be mapped to obtain effective boost terminal states. The same reachable area may correspond to multiple sets of boost terminal states. However, the reachable area boundary parameters caused by the deviation of the boost terminal states are unique. Therefore, it is possible to construct the mapping relationship to realize the transformation from the reachable area boundary to the terminal state increments of the boost phase. Thus, the networks of altitude, velocity, and flight path angle as shown in Figure7 are established.Figure 7
Relational mapping diagram of deep neural networks.Generally, the neural networks with 1~2 hidden layers can fit the nonlinear smooth mapping with arbitrary accuracy, and the ability of neural networks to describe complex learning can be increased by introducing the third hidden layer [25]. For the training of neural networks with three parameters of the reachable area boundary as input and one terminal state increment of the boost phase as output, the number of hidden layers was set to 3 in this study. For the hyperparameter setting of other hidden layer neural networks, some meaningful manual optimization rules were given in Nielsen’s works [26]. In literature [27, 28], an automatic optimization method of hyperparameters using grid search technology was introduced. Therefore, this paper only briefly gives the process of determining the number of hidden layers: firstly, the number of nodes in each hidden layer is set to 2~20; then, the number of hidden layer nodes is increased in turn for network training; finally, the network structure with the least total error on the verification is recorded as the optimal neural network structure.


### 4.3. Energy Equivalence Principles

It is necessary to adopt the optimal adjustment strategy of the reachable area according to the influence of the boost terminal states considering limited fuel. Under fixed time constraints, the trajectory optimization was carried out with the boost terminal point’s altitude, velocity, and flight path angle as performance indexes. By taking the obtained boost terminal states as the initial condition for the glider, the reachable area boundaries under different energy consumption conditions (100 kg, 200 kg, and 300 kg fuel increase, respectively) as shown in Figure8 are obtained by the method in Section 3.1.Figure 8
Effect of terminal states of the boost phase on the reachable area.Under the same energy consumption conditions, the reachable area obtained by raising the velocity of the boost terminal was larger than the reachable area obtained by raising the flight path angle of the boost terminal, and both were larger than the reachable area obtained by increasing the height of the boost terminal. Therefore, on the basis of getting the new reachable area boundary parameters, the terminal states of the boost phase can be obtained through DNNs described in Section4.2. The selection priority of the mapping networks is velocity networks>flight path angle networks>altitude networks.Although the selection priority can maximize the reachable area, the terminal states of the boost phase must be within a certain limit due to the debris drop point constraints, separation conditions, and the limitation of the glider’s navigation communication window. Assume that the separation window constraints of the booster based on the nominal terminal states areΔhmax,Δvmax,Δγmax. When the boost terminal states obtained by DNN exceed the separation window constraints, the energy equivalence principles are adopted to obtain the new boost terminal states.The method described in Section3 can be used to obtain the relationship between the boost terminal states and the reachable area boundary parameters, as shown in Table 4.Table 4
Relationship between boost terminal states and reachable area boundary parameters.
StatesStates changeVariation of boundary curve parametersHeight (km)0Δa=1.0068719×10−9,Δb=4.7115307×10−8,Δymax=010.4Δa=8.3883981×10−3,Δb=1.0403449×10−2,Δymax=3.3730000×10−317.3Δa=1.0925126×10−2,Δb=2.0196111×10−2,Δymax=2.6560000×10−323.25Δa=5.5262781×10−3,Δb=2.7402044×10−2,Δymax=4.5460000×10−3Velocity (m/s)0Δa=1.0068719×10−9,Δb=4.7115307×10−8,Δymax=069.2Δa=4.2817309×10−2,Δb=4.9985394×10−2,Δymax=1.8905000×10−2142Δa=−9.6140915×10−2,Δb=1.0372763×10−1,Δymax=4.4436000×10−2216.3Δa=−1.4838703×10−1,Δb=1.5801915×10−1,Δymax=6.5828000×10−2Flight path angle (°)4.25Δa=8.9585021×10−2,Δb=6.2255655×10−2,Δymax=−3.6950000×10−36.5Δa=−1.5504815×10−1,Δb=8.8706775×10−2,Δymax=−9.7200000×10−37.8Δa=2.1084387×10−1,Δb=1.0807930×10−1,Δymax=−1.0417000×10−28.9Δa=2.1735414×10−1,Δb=1.0080363×10−1,Δymax=−2.1716000×10−2The contents of the energy equivalence principles are as follows. ① After the reachable area boundary parameters are obtained by the constraint transformation rules, the velocity networks can be used to obtain the boost terminal velocity incrementΔv. If Δv≤Δvmax, the final boost terminal state constraints are h0,v0+Δv,γ0; if Δv>Δvmax, proceed to the next step. ② Let δΔv=Δv−Δvmax and interpolate Table 4 with δΔv to obtain the reachable area boundary parameter increment Δa,b,ymaxΔv; then, the flight path angle neural networks can be used to obtain the new boost terminal increment of the flight path angle Δγ. If Δγ≤Δγmax, the final boost terminal state constraints are h0,v0+Δvmax,γ0+Δγ; if Δγ>Δγmax, proceed to the next step. ③ Let δΔγ=Δγ−Δγmax and interpolate Table 4 with δΔγ to obtain the increment of reachable area boundary parameters Δa,b,ymaxΔγ; then, the height neural networks are used to obtain the new boost terminal height increment Δh. If Δh≤Δhmax, the final boost terminal state constraints are h0+Δh,v0+Δvmax,γ0+Δγmax; if Δh>Δhmax, the constraint window is exceeded, and the boost terminal states are set to the maximum separation window constraints h0+Δhmax,v0+Δvmax,γ0+Δγmax.


## 4.1. Constraint Transformation Rules

This section describes the method to transform NRACs to generate new coverage envelopes about the reachable area. As shown in Figure5, the projection of the boost terminal point on the ground is the origin O, and the coordinate system is established in the horizontal plane, with the axis Ox in the same direction as the projection of the initial velocity direction v0 on the ground. The original terminal states of the boost phase are h0,v0,γ0; x=a0y2+b0−ymax≤y≤ymax is the reachable area boundary curve generated by h0,v0,γ0; the nonregular reachable area is Q; the distance from the farthest point in Q to the coordinate origin is R; the angle between lL and lR formed by the left and right boundary points of Q and the coordinate origin O is Δξ; the centerline of Δξ is the new axis Ox′; the angle between Ox′ and Ox is Δψt; LT:x=ay2+b denotes the new reachable area boundary curve, so the key of the issue is to solve the parameters of LT.Figure 5
Description of the coverage requirements of the nonregular reachable area.Whether the original terminal states of the boost phase are appropriate can be judged by the relative position relationship between the boundary points ofQ and the boundary of the original reachable area.
(9)Judgment1:a0RsinΔψt2+b0>RcosΔψt,Judgment2:a0l0sinΔψt−12Δξ+Δξ0n2+b0>l0cosΔψt−12Δξ+Δξ0n,Judgment3:a0l1sinΔψt−12Δξ+Δξ1n2+b0>l0cosΔψt−12Δξ+Δξ1n⋯Judgmenti:a0li−2sinΔψt−12Δξ+Δξi−2n2+b0>li−2cosΔψt−12Δξ+Δξi−2n,where n is the number of divisions when judging the boundary for Q; li is the distance from the coordinate origin to the i‐th boundary judgment point. The n+2 judgment conditions must be satisfied simultaneously to determine whether the original terminal states of the boost phase are appropriate. The number of divisions can be obtained according to the description accuracy of the reachable area and the shape of Q.When the original terminal states of the boost phase do not meet NRACs, the change of initial heading angle can bring great convenience to the adjustment of the reachable area. However, the reachable area adjusted by the initial flight direction may still not meet the coverage requirements. Therefore, the judgments need to be reconducted in the newOx′y′-coordinate system.
(10)Judgment1:b0>R,Judgment2:a0l0sin−12Δξ+Δξ0n2+b0>l0cos−12Δξ+Δξ0n,Judgment3:a0l1sin−12Δξ+Δξ1n2+b0>l0cos−12Δξ+Δξ1n,⋯Judgmenti:a0li−2sin−12Δξ+Δξi−2n2+b0>li−2cos−12Δξ+Δξi−2n.When then+2 judgment conditions cannot be satisfied simultaneously, the new boundary parameters can be adjusted and determined according to the range of values of the boundary parameters.
(11)amin≤a≤amax,bmin≤b≤bmax,ymaxmin≤ymax≤ymaxmax.As a supporting parameter, the parameterymax varies inversely and equivalently with the adjustment of a. The adjustment strategies for the three parameters are shown in Figure 6.Figure 6
Adjustment strategies for the reachable area boundary parameters.


## 4.2. Relational Mapping of Deep Neural Networks

The relationship between the reachable area boundary parameters and the boost terminal states needs to be mapped to obtain effective boost terminal states. The same reachable area may correspond to multiple sets of boost terminal states. However, the reachable area boundary parameters caused by the deviation of the boost terminal states are unique. Therefore, it is possible to construct the mapping relationship to realize the transformation from the reachable area boundary to the terminal state increments of the boost phase. Thus, the networks of altitude, velocity, and flight path angle as shown in Figure7 are established.Figure 7
Relational mapping diagram of deep neural networks.Generally, the neural networks with 1~2 hidden layers can fit the nonlinear smooth mapping with arbitrary accuracy, and the ability of neural networks to describe complex learning can be increased by introducing the third hidden layer [25]. For the training of neural networks with three parameters of the reachable area boundary as input and one terminal state increment of the boost phase as output, the number of hidden layers was set to 3 in this study. For the hyperparameter setting of other hidden layer neural networks, some meaningful manual optimization rules were given in Nielsen’s works [26]. In literature [27, 28], an automatic optimization method of hyperparameters using grid search technology was introduced. Therefore, this paper only briefly gives the process of determining the number of hidden layers: firstly, the number of nodes in each hidden layer is set to 2~20; then, the number of hidden layer nodes is increased in turn for network training; finally, the network structure with the least total error on the verification is recorded as the optimal neural network structure.


## 4.3. Energy Equivalence Principles

It is necessary to adopt the optimal adjustment strategy of the reachable area according to the influence of the boost terminal states considering limited fuel. Under fixed time constraints, the trajectory optimization was carried out with the boost terminal point’s altitude, velocity, and flight path angle as performance indexes. By taking the obtained boost terminal states as the initial condition for the glider, the reachable area boundaries under different energy consumption conditions (100 kg, 200 kg, and 300 kg fuel increase, respectively) as shown in Figure8 are obtained by the method in Section 3.1.Figure 8
Effect of terminal states of the boost phase on the reachable area.Under the same energy consumption conditions, the reachable area obtained by raising the velocity of the boost terminal was larger than the reachable area obtained by raising the flight path angle of the boost terminal, and both were larger than the reachable area obtained by increasing the height of the boost terminal. Therefore, on the basis of getting the new reachable area boundary parameters, the terminal states of the boost phase can be obtained through DNNs described in Section4.2. The selection priority of the mapping networks is velocity networks>flight path angle networks>altitude networks.Although the selection priority can maximize the reachable area, the terminal states of the boost phase must be within a certain limit due to the debris drop point constraints, separation conditions, and the limitation of the glider’s navigation communication window. Assume that the separation window constraints of the booster based on the nominal terminal states areΔhmax,Δvmax,Δγmax. When the boost terminal states obtained by DNN exceed the separation window constraints, the energy equivalence principles are adopted to obtain the new boost terminal states.The method described in Section3 can be used to obtain the relationship between the boost terminal states and the reachable area boundary parameters, as shown in Table 4.Table 4
Relationship between boost terminal states and reachable area boundary parameters.
StatesStates changeVariation of boundary curve parametersHeight (km)0Δa=1.0068719×10−9,Δb=4.7115307×10−8,Δymax=010.4Δa=8.3883981×10−3,Δb=1.0403449×10−2,Δymax=3.3730000×10−317.3Δa=1.0925126×10−2,Δb=2.0196111×10−2,Δymax=2.6560000×10−323.25Δa=5.5262781×10−3,Δb=2.7402044×10−2,Δymax=4.5460000×10−3Velocity (m/s)0Δa=1.0068719×10−9,Δb=4.7115307×10−8,Δymax=069.2Δa=4.2817309×10−2,Δb=4.9985394×10−2,Δymax=1.8905000×10−2142Δa=−9.6140915×10−2,Δb=1.0372763×10−1,Δymax=4.4436000×10−2216.3Δa=−1.4838703×10−1,Δb=1.5801915×10−1,Δymax=6.5828000×10−2Flight path angle (°)4.25Δa=8.9585021×10−2,Δb=6.2255655×10−2,Δymax=−3.6950000×10−36.5Δa=−1.5504815×10−1,Δb=8.8706775×10−2,Δymax=−9.7200000×10−37.8Δa=2.1084387×10−1,Δb=1.0807930×10−1,Δymax=−1.0417000×10−28.9Δa=2.1735414×10−1,Δb=1.0080363×10−1,Δymax=−2.1716000×10−2The contents of the energy equivalence principles are as follows. ① After the reachable area boundary parameters are obtained by the constraint transformation rules, the velocity networks can be used to obtain the boost terminal velocity incrementΔv. If Δv≤Δvmax, the final boost terminal state constraints are h0,v0+Δv,γ0; if Δv>Δvmax, proceed to the next step. ② Let δΔv=Δv−Δvmax and interpolate Table 4 with δΔv to obtain the reachable area boundary parameter increment Δa,b,ymaxΔv; then, the flight path angle neural networks can be used to obtain the new boost terminal increment of the flight path angle Δγ. If Δγ≤Δγmax, the final boost terminal state constraints are h0,v0+Δvmax,γ0+Δγ; if Δγ>Δγmax, proceed to the next step. ③ Let δΔγ=Δγ−Δγmax and interpolate Table 4 with δΔγ to obtain the increment of reachable area boundary parameters Δa,b,ymaxΔγ; then, the height neural networks are used to obtain the new boost terminal height increment Δh. If Δh≤Δhmax, the final boost terminal state constraints are h0+Δh,v0+Δvmax,γ0+Δγmax; if Δh>Δhmax, the constraint window is exceeded, and the boost terminal states are set to the maximum separation window constraints h0+Δhmax,v0+Δvmax,γ0+Δγmax.


## 5. Boost-Phase Trajectory Planning

After obtaining the terminal state constraints of the boost phase, the trajectory from the initial point to the new terminal point needs to be replanned. Unlike the trajectory tracking with a known reference trajectory [29], the boost-phase trajectory is replanned under the condition that the initial value is not ideal or there is no reference trajectory in this study. At the same time, the trajectory planning method needs to overcome the impacts of initial deviation and environmental disturbance conditions and needs to meet the process constraints such as states and control and the given terminal states constraints.The Gaussian pseudospectral method can be used to realize fast and accurate trajectory replanning to meet the requirements of low sensitivity to initial values, good convergence, and high accuracy. Aiming at the “Dimensional Catastrophe” of Jacobian and Hessian matrices and the slow convergence speed when solving nonsmooth problems, the hp-type adaptive finite element method (hp-FEM) is used, which allows the segment length and the order of basis function to change adaptively on the mesh profile at the same time.By taking the longest range of the boost phase is taken as the performance index, the boost-phase trajectory planning problem is converted into an optimal control problem:(12)J=−Re⋅arccoscosφfcosφ0cosθf−θ0+sinφfsinφ0,ẋt=fxt,ut,t,t∈t0,tf,ϕxt0,t0,xtf,tf=0,Cxt,ut,t≤0,t∈t0,tf,where xt denotes the states; ut denotes the control variables; ẋt=fxt,ut,t denotes the dynamic equations corresponding to Equation (1); t0,tf is the time interval; ϕxt0,t0,xtf,tf denotes the boundary constraints; Cxt,ut,t denotes the path constraints; θ0,φ0 is the latitude and longitude of the launch point; θf,φf is the latitude and longitude of the terminal location of the boost phase.The control variables are the attack angle and sideslip angle:(13)ut=α,βBoundary constraints are as follows:(14)hxtf,tf−hf=0,vxtf,tf−vf=0,γxtf,tf−γf=0,ψxtf,tf−ψf=0,where hf, vf, γf, and ψf denote the terminal altitude, velocity, flight path angle, and heading angle, respectively.Path constraints are as follows:(15)xmin≤ϕxt,t≤xmax,αmin≤α≤αmax,α̇≤α̇max,βmin≤β≤βmax,β̇≤β̇max.The time, states, and control variables must be transformed and discretized.Xτ denotes the discretized states, and Uτ denotes the control variables. The time interval t0,tf is converted into [-1,1].
(16)τ=2ttf−t0−tf+t0tf−t0.TakingK‐th Legendre-Gauss points τ1,⋯,τK and τ0=−1 as nodes, the K+1‐th Lagrange interpolation polynomials are formed as basis functions to approximate the states.
(17)xτ≈Xτ=∑i=0KLiτxτiLiτ=∏j=0,j≠iKτ−τjτi−τj.,The approximation of control variables are as follows:(18)uτ≈Uτ=∑i=1KL~iτuτi,L~iτ=∏j=1,j≠iKτ−τjτi−τj.Furthermore, the kinetic differential equations can be transformed into algebraic constraints.(19)ẋτk≈Ẋτk=∑i=0KL̇iτkXτi=∑i=0KDkiτkXτi,∑i=0KDkiτkXτi−tf−t02fXτk,Uτk,τk;t0,tf=0,k=1,⋯,K,whereDki is the differential matrix derived from the Legendre polynomial; PKτ is the root of the K‐th Legendre polynomial [20].
(20)Dki=L̇iτk=1+τkṖKτk+PKτkτk−τi1+τiṖKτi+PKτi,i≠k,1+τiP¨Kτi+2ṖKτi21+τiṖKτi+PKτi,i=k,PKτ=12KK!dKdτKτ2−1K.Discrete approximation of path constraints are as follows:(21)xmin≤Xτ≤xmax,αmin;βmin≤Uτ≤αmax;βmax,U̇τj=∑i=1KL~̇iτjUτi≤α̇max;β̇max.Discrete approximation of boundary constraints are as follows:(22)hXτf,τf−hf=0,vXτf,τf−vf=0,γXτf,τf−γf=0,ψXτf,τf−ψf=0.Since the terminal states are not included in the state approximation, a discrete approximation can be performed using Gaussian integration to satisfy the dynamics constraints.(23)Xτf−Xτ0−tf−t02∑k=1KwkfXτk,Uτk,τk;t0,tf=0,wk=∫−11L~kτdτ=21−τk2ṖKτk2.The optimal control problem has been transformed into a discrete nonlinear programming problem, which can be solved using the large-scale sequential quadratic programming-based SNOPT solver.There have been a great number of research achievements on using the hp-adaptive pseudospectral method to solve the trajectory planning problem [30–32]. The convergence of the Gaussian pseudospectral method has been studied and described in reference [33]. Therefore, this paper will not repeat this part but only discuss the convergence range of the deviation conditions in the simulation section.


## 6. Simulation and Analysis

### 6.1. Simulation Settings

The simulation configuration is Intel (R) Core (TM) i7-10510U. The general parameters of the vehicle, mission parameters, disturbances, and deviation conditions required for the simulation are given in Tables5–7.Table 5
General parameters of the vehicle.
ItemValue (unit)Boost phaseInitial mass10000 kgConsumption per second60 kg/sSpecific impulse300 sPneumatic reference area1 m2Fuel mass5850 kgGliding phaseMass1000 kgPneumatic reference area0.5 m2Table 6
Mission parameters.
ItemValue (unit)Air-based launch pointθ=0°,φ=0°The initial velocity of the boost phaseV=3101m/s,θv=9.31°The initial position of the boost phasex=635000m,y=−3000m,z=10000mNominal terminal states of the boost phaseh=70km,v=5000m/s,γ=−1°,ψ=111°Control constraints−20°≤α≤20°, −10°≤β≤10°Control rate constraintsα̇≤10°/s, β̇≤5°/sSeparation window constraintsΔvmax=100m/s,Δγmax=3°,Δhmax=10kmTerminal states of the gliding phasehf=20km,vf=1000m/sRequired coverage areaP1:54°,18°,P2:59°,18°,P3:56°,13°,P4:52°,11°Table 7
Disturbances and deviation conditions.
ItemDeviation valueMaximum energyMinimum energyInitial state deviationsPosition (m)[2000, 2000, 1000][-2000, -2000, -1000]Velocity (m/s)[20, 20, 10][-20, -20, -10]Power system deviationFuel mass (kg)50-50Environmental disturbancesResistance coefficient (%)-1515Lift coefficient (%)15-15Atmospheric density (%)-1010


### 6.2. Simulation Results

(1)
Training results of deep neural networksThree different DNNs were trained using the parameterized reachable area boundary library in Section3.2 as the samples. From the single hidden layer network structure of three neurons, the depth and breadth of neural network are gradually increased until the accuracy requirements are met. The training, verification, and test results are shown in Figures 9–11.Figure 9
Training results of velocity neural networks.Figure 10
Training results of flight path angle neural networks.Figure 11
Training results of high neural networks.The simulation results show that all three deep neural networks have achieved convergence, with training and testing accuracy higher than 99%. The generated demand reachable area boundary parameter increments areΔa=0,Δb=0.055,Δymax=0 based on the previous geometric relations.
(2)
Simulation of boost-phase trajectory planningTo verify the effectiveness of the proposed trajectory planning method and analyze the algorithm’s performance, NSGA II as a widely used evolutionary algorithm was used in simulation experiments for comparison. The principle of NSGA II is referred to in [34].Optimization objectives(24)JN1=Δhf,JN2=Δvf,JN3=Δγf,JN4=Δψf,JN5=−Lrange,where Δhf, Δvf, Δγf, and Δψf are the terminal states deviations; Lrange denotes the terminal range.Control variables(25)uNi=αNi,i=1,2,⋯,5,uN6=β̇N,where αNi denotes i‐th of the attack angle at linear interpolation node; β̇N denotes the change rate of sideslip angle.Constraints(26)−20°≤uNi≤20°,i=1,2,⋯,5,−5°/s≤uNi≤5°/s,i=6.Other main parameters are shown in Table8.Table 8
Parameters of NSGA II.
ItemValuePopulationSize200Generations80StallGenLimit40CrossoverFraction0.8MigrationFraction0.4ParetoFraction0.3The terminal constraints of the boost phase obtained by using neural networks areΔv=100m/s, Δγ=0.8°, Δh=0km, Δψ=−21.3°. “hpPM” represents the planning result of the hp-adaptive pseudospectral method, “verification” denotes the integral verification result of the hp-adaptive pseudospectral method, and “ NSGA II” represents the planning result of the improved nondominated sorting genetic algorithm. The trajectory planning curves are shown in Figures 12–21, and the comparison results of two methods are shown in Table 9.Figure 12
Height.Figure 13
Velocity.Figure 14
Flight path angle.Figure 15
Heading angle.Figure 16
Attack angle.Figure 17
Sideslip angle.Figure 18
Change rate of attack angle.Figure 19
Change rate of sideslip angle.Figure 20
Longitudinal range.Figure 21
Pareto front.Table 9
Comparison results.
ItemhpPMNSGA IITerminal deviationsHeight (m)0.076.0Velocity (m/s)0.09.6Flight path angle (°)0.0−3.4×10−2Heading angle (°)1.0×10−3−5.0×10−3Performance indexMaximum range (km)988.5987.1Path constraintsMaximum attack angle (°)20.020.0Maximum change rate of attack angle (°/s)5.834.2Maximum sideslip angle (°)10.07.3Maximum change rate of sideslip angle (°/s)1.43.0×10−2TimelinessCalculation time (s)2.8318.3The results show that ① the boost-phase flight trajectory obtained by the hpPM that met the requirements of dynamic integration verification; ② the terminal planning accuracy of hpPM was higher than that of NSGA II; ③ the trajectory obtained by hpPM had a larger range than that obtained by NSGA II; ④ the change rate of control angle planned by hpPM was more stable than that by NSGA II in longitudinal profile; the control angle and the change rate of control angle were the same in transverse profile; ⑤ the timeliness of hpPM was significantly better than that of NSGA II.(3)
Influence of disturbances and deviationsDue to the influence of NRACs, the boost-phase trajectory planning with fixed working time has been a challenge. To analyze the impacts of initial state deviations, power system deviation, and environmental disturbances on hpPM, it is necessary to conduct simulation experiments under the limit energy deviation as shown in Table7. The trajectory planning curves are shown in Figures 22–31, and the results of model uncertainty and disturbances are shown in Table 10.Figure 22
Height.Figure 23
Velocity.Figure 24
Flight path angle.Figure 25
Heading angle.Figure 26
Attack angle.Figure 27
Sideslip angle.Figure 28
Change rate of attack angle.Figure 29
Change rate of sideslip angle.Figure 30
Longitudinal range.Figure 31
Mass.Table 10
Results of model uncertainty and disturbances.
ItemMaximum energyMinimum energyTerminal deviationHeight (m)500.00.0Velocity (m/s)0.0-71.0Flight path angle (°)0.00.0Heading angle (°)0.00.0Performance indexMaximum range (km)996.0978.5Path constraintsMaximum attack angle (°)20.020.0Maximum change rate of attack angle (°/s)10.010.0Maximum sideslip angle (°)10.010.0Maximum change rate of sideslip angle (°/s)5.02.8TimelinessCalculation time (s)8.614.9Due to limited fuel, the adjustment ability of the boost-glide vehicle is limited. Therefore, a limited energy deviation of trajectory planning failure is of great importance. The minimum energy deviation will seriously reduce the flight maneuverability, making the boost-glide vehicle barely reach the designated mission point with the loss of working performance; the maximum energy deviation will make the boost-glide vehicle produce excess flight maneuverability, leading to a large oscillation of attitude angle and instability of flight.The results show that ① for the maximum energy deviation, the velocity, flight path angle, and heading angle of boost-glide vehicle met the terminal indexes, but the altitude and range were greater than the nominal case, and an over large change rate of control angle was produced. If the magnitude of energy deviations continued to increase, the boost-glide vehicle would lose stability easily or exceed the separation window constraint; ② for the minimum energy deviation, the altitude, flight path angle, and heading angle of the boost-glide vehicle fully met the terminal indexes, but the velocity and range were less than the nominal case. If the magnitude of energy deviations continued to increase, the reachable area of the boost-glide vehicle would be seriously reduced. Therefore, the disturbances and deviations listed in Table7 can be considered as the failure range of the trajectory planned by hpPM under the mission specified in Table 6.
(4)
The verification simulation of the reachable areaAfter obtaining the deviations of the range, height, and velocity of the boost terminal point, the distribution of reachable area can be calculated, and the effectiveness of the proposed trajectory planning method considering NRACs can be further analyzed. In the comparison of reachable area distribution, “nominal reachable area” refers to the maneuvering range of the boost-glider vehicle without any adjustments; “demand boundary” refers to the maneuvering range meeting NRACs obtained by the constraint transformation rules; “maximum/minimum energy” represents the maneuvering range under the influence of disturbances and deviations. The results are shown in Figure32.Figure 32
Comparison of the reachable area.The results show that ① the nominal reachable area could not meet the requirements of NRACs, but the demand boundary transformed by the method proposed in this paper completely covered the area specified by NRACs; ② the actual maneuvering range of the boost-glide vehicle still met NRACs under the influence of disturbances and deviations; ③ the average deviation between the actual reachable area boundary and the demand boundary under the limit energy deviation was 3.9%. The above results indicated that the proposed trajectory planning method considering NRACs was sufficiently capable of adjusting the reachable area and had high trajectory planning accuracy and good deviation adaptability.


## 6.1. Simulation Settings

The simulation configuration is Intel (R) Core (TM) i7-10510U. The general parameters of the vehicle, mission parameters, disturbances, and deviation conditions required for the simulation are given in Tables5–7.Table 5
General parameters of the vehicle.
ItemValue (unit)Boost phaseInitial mass10000 kgConsumption per second60 kg/sSpecific impulse300 sPneumatic reference area1 m2Fuel mass5850 kgGliding phaseMass1000 kgPneumatic reference area0.5 m2Table 6
Mission parameters.
ItemValue (unit)Air-based launch pointθ=0°,φ=0°The initial velocity of the boost phaseV=3101m/s,θv=9.31°The initial position of the boost phasex=635000m,y=−3000m,z=10000mNominal terminal states of the boost phaseh=70km,v=5000m/s,γ=−1°,ψ=111°Control constraints−20°≤α≤20°, −10°≤β≤10°Control rate constraintsα̇≤10°/s, β̇≤5°/sSeparation window constraintsΔvmax=100m/s,Δγmax=3°,Δhmax=10kmTerminal states of the gliding phasehf=20km,vf=1000m/sRequired coverage areaP1:54°,18°,P2:59°,18°,P3:56°,13°,P4:52°,11°Table 7
Disturbances and deviation conditions.
ItemDeviation valueMaximum energyMinimum energyInitial state deviationsPosition (m)[2000, 2000, 1000][-2000, -2000, -1000]Velocity (m/s)[20, 20, 10][-20, -20, -10]Power system deviationFuel mass (kg)50-50Environmental disturbancesResistance coefficient (%)-1515Lift coefficient (%)15-15Atmospheric density (%)-1010


## 6.2. Simulation Results

(1)
Training results of deep neural networksThree different DNNs were trained using the parameterized reachable area boundary library in Section3.2 as the samples. From the single hidden layer network structure of three neurons, the depth and breadth of neural network are gradually increased until the accuracy requirements are met. The training, verification, and test results are shown in Figures 9–11.Figure 9
Training results of velocity neural networks.Figure 10
Training results of flight path angle neural networks.Figure 11
Training results of high neural networks.The simulation results show that all three deep neural networks have achieved convergence, with training and testing accuracy higher than 99%. The generated demand reachable area boundary parameter increments areΔa=0,Δb=0.055,Δymax=0 based on the previous geometric relations.
(2)
Simulation of boost-phase trajectory planningTo verify the effectiveness of the proposed trajectory planning method and analyze the algorithm’s performance, NSGA II as a widely used evolutionary algorithm was used in simulation experiments for comparison. The principle of NSGA II is referred to in [34].Optimization objectives(24)JN1=Δhf,JN2=Δvf,JN3=Δγf,JN4=Δψf,JN5=−Lrange,where Δhf, Δvf, Δγf, and Δψf are the terminal states deviations; Lrange denotes the terminal range.Control variables(25)uNi=αNi,i=1,2,⋯,5,uN6=β̇N,where αNi denotes i‐th of the attack angle at linear interpolation node; β̇N denotes the change rate of sideslip angle.Constraints(26)−20°≤uNi≤20°,i=1,2,⋯,5,−5°/s≤uNi≤5°/s,i=6.Other main parameters are shown in Table8.Table 8
Parameters of NSGA II.
ItemValuePopulationSize200Generations80StallGenLimit40CrossoverFraction0.8MigrationFraction0.4ParetoFraction0.3The terminal constraints of the boost phase obtained by using neural networks areΔv=100m/s, Δγ=0.8°, Δh=0km, Δψ=−21.3°. “hpPM” represents the planning result of the hp-adaptive pseudospectral method, “verification” denotes the integral verification result of the hp-adaptive pseudospectral method, and “ NSGA II” represents the planning result of the improved nondominated sorting genetic algorithm. The trajectory planning curves are shown in Figures 12–21, and the comparison results of two methods are shown in Table 9.Figure 12
Height.Figure 13
Velocity.Figure 14
Flight path angle.Figure 15
Heading angle.Figure 16
Attack angle.Figure 17
Sideslip angle.Figure 18
Change rate of attack angle.Figure 19
Change rate of sideslip angle.Figure 20
Longitudinal range.Figure 21
Pareto front.Table 9
Comparison results.
ItemhpPMNSGA IITerminal deviationsHeight (m)0.076.0Velocity (m/s)0.09.6Flight path angle (°)0.0−3.4×10−2Heading angle (°)1.0×10−3−5.0×10−3Performance indexMaximum range (km)988.5987.1Path constraintsMaximum attack angle (°)20.020.0Maximum change rate of attack angle (°/s)5.834.2Maximum sideslip angle (°)10.07.3Maximum change rate of sideslip angle (°/s)1.43.0×10−2TimelinessCalculation time (s)2.8318.3The results show that ① the boost-phase flight trajectory obtained by the hpPM that met the requirements of dynamic integration verification; ② the terminal planning accuracy of hpPM was higher than that of NSGA II; ③ the trajectory obtained by hpPM had a larger range than that obtained by NSGA II; ④ the change rate of control angle planned by hpPM was more stable than that by NSGA II in longitudinal profile; the control angle and the change rate of control angle were the same in transverse profile; ⑤ the timeliness of hpPM was significantly better than that of NSGA II.(3)
Influence of disturbances and deviationsDue to the influence of NRACs, the boost-phase trajectory planning with fixed working time has been a challenge. To analyze the impacts of initial state deviations, power system deviation, and environmental disturbances on hpPM, it is necessary to conduct simulation experiments under the limit energy deviation as shown in Table7. The trajectory planning curves are shown in Figures 22–31, and the results of model uncertainty and disturbances are shown in Table 10.Figure 22
Height.Figure 23
Velocity.Figure 24
Flight path angle.Figure 25
Heading angle.Figure 26
Attack angle.Figure 27
Sideslip angle.Figure 28
Change rate of attack angle.Figure 29
Change rate of sideslip angle.Figure 30
Longitudinal range.Figure 31
Mass.Table 10
Results of model uncertainty and disturbances.
ItemMaximum energyMinimum energyTerminal deviationHeight (m)500.00.0Velocity (m/s)0.0-71.0Flight path angle (°)0.00.0Heading angle (°)0.00.0Performance indexMaximum range (km)996.0978.5Path constraintsMaximum attack angle (°)20.020.0Maximum change rate of attack angle (°/s)10.010.0Maximum sideslip angle (°)10.010.0Maximum change rate of sideslip angle (°/s)5.02.8TimelinessCalculation time (s)8.614.9Due to limited fuel, the adjustment ability of the boost-glide vehicle is limited. Therefore, a limited energy deviation of trajectory planning failure is of great importance. The minimum energy deviation will seriously reduce the flight maneuverability, making the boost-glide vehicle barely reach the designated mission point with the loss of working performance; the maximum energy deviation will make the boost-glide vehicle produce excess flight maneuverability, leading to a large oscillation of attitude angle and instability of flight.The results show that ① for the maximum energy deviation, the velocity, flight path angle, and heading angle of boost-glide vehicle met the terminal indexes, but the altitude and range were greater than the nominal case, and an over large change rate of control angle was produced. If the magnitude of energy deviations continued to increase, the boost-glide vehicle would lose stability easily or exceed the separation window constraint; ② for the minimum energy deviation, the altitude, flight path angle, and heading angle of the boost-glide vehicle fully met the terminal indexes, but the velocity and range were less than the nominal case. If the magnitude of energy deviations continued to increase, the reachable area of the boost-glide vehicle would be seriously reduced. Therefore, the disturbances and deviations listed in Table7 can be considered as the failure range of the trajectory planned by hpPM under the mission specified in Table 6.
(4)
The verification simulation of the reachable areaAfter obtaining the deviations of the range, height, and velocity of the boost terminal point, the distribution of reachable area can be calculated, and the effectiveness of the proposed trajectory planning method considering NRACs can be further analyzed. In the comparison of reachable area distribution, “nominal reachable area” refers to the maneuvering range of the boost-glider vehicle without any adjustments; “demand boundary” refers to the maneuvering range meeting NRACs obtained by the constraint transformation rules; “maximum/minimum energy” represents the maneuvering range under the influence of disturbances and deviations. The results are shown in Figure32.Figure 32
Comparison of the reachable area.The results show that ① the nominal reachable area could not meet the requirements of NRACs, but the demand boundary transformed by the method proposed in this paper completely covered the area specified by NRACs; ② the actual maneuvering range of the boost-glide vehicle still met NRACs under the influence of disturbances and deviations; ③ the average deviation between the actual reachable area boundary and the demand boundary under the limit energy deviation was 3.9%. The above results indicated that the proposed trajectory planning method considering NRACs was sufficiently capable of adjusting the reachable area and had high trajectory planning accuracy and good deviation adaptability.


## 7. Conclusion

Aiming at the boost-phase trajectory planning under NRACs, this paper studies the parametric description of the reachable area, the reachable area constraint transformation, and the boost-phase trajectory planning. The NRACs were successfully transformed into the terminal state constraints of the boost phase, and the hpPM was used to complete the trajectory planning. The main conclusions are as follows.(1)
The reachable area boundary of the boost-glide vehicle can be described and represented in a simple parametric form. The quadratic form is a feasible solution considering the accuracy and parameter complexity(2)
There is a clear relationship between the reachable area boundary described by parameterization and the boost terminal state change. The DNNs can be used for model learning and parameter mapping, with accuracy of model training and testing greater than 99%. In addition, the changes of boost terminal states pose a significant influence on the reachable area, that is, the influence of velocity was greater than that of flight path angle, and the influence of altitude is the most insignificant(3)
Using the constraint transformation rules, it is possible to quickly determine the reachable area boundary parameters that meet NRACs after changing the initial flight direction of the boost-glide vehicle through the air-based launch mode or the promotion of the front stage(4)
Based on the deep neural networks model and the energy equivalence principles, the appropriate boost terminal states under the separation window constraints can be obtained(5)
The boost-phase trajectory planning method based on the hpPM has good timeliness, high planning accuracy, and strong deviation adaptabilitySince this paper focuses on the boost-phase trajectory planning method considering NRACs, the determination of the reachable area of the boost-glide vehicle is mainly depended on the fast calculation method described in Section3.1 in spite of the influences of the deviations and disturbances of the gliding phase. Future research will focus on the influences of combined deviations in different stages.



---

*Source: 9443050-2022-06-26.xml*