Aircraft trajectory calculations include dynamic and kinematic models. The equations describing the motion parameters of the aircraft centroid include: { R ˙ = V ⁡ sin ⁡ γ θ ˙ = V ⁡ cos ⁡ γ ⁡ sin ⁡ ψ V R ⁡ cos ⁡ ϕ ϕ ˙ = V ⁡ cos ⁡ γ ⁡ cos ⁡ ψ V R V ˙ = P ⁡ cos ⁡ α ⁡ cos ⁡ β − D m − g ⁡ sin ⁡ γ γ ˙ = 1 m V [ P ( sin ⁡ α ⁡ cos ⁡ γ V + cos ⁡ α ⁡ cos ⁡ β ⁡ sin ⁡ γ V ) + L ⁡ cos ⁡ γ V − Z ⁡ sin ⁡ γ V ] − g V cos ⁡ γ ψ V = − 1 m V ⁡ cos ⁡ γ [ P ( sin ⁡ α ⁡ sin ⁡ γ V − cos ⁡ α ⁡ sin ⁡ β ⁡ cos ⁡ γ V ) + L ⁡ sin ⁡ γ V + Z ⁡ cos ⁡ γ V ] , (1) where R is the distance from the aircraft to the earth’s center, V is the aircraft speed, θ and ϕ are the longitude and latitude of the aircraft, respectively, γ is the trajectory inclination, ψ V is the trajectory deflection angle, γ V is the speed inclination angle, α and β are the attack angle and sideslip angle of the aircraft respectively, L, D, Z, and P are lift, drag, lateral force, and engine thrust, respectively.

The aerodynamic force acting on the aircraft is the functional relationship of flight speed Ma, height H, attitude angles α , β , etc., and control surface deflection angles δ x , δ y , δ z , etc., which can be expressed as: C L , D , Z = f ( M a , H , α , β , δ x , δ y , δ z ) , (2)

In addition, the trajectory calculation also needs the geometric relationship between the angles in formula (1), atmospheric model, a control system loop model, etc.

For hypersonic vehicles, they go through different flight stages, from ground zero speed take-off to high-altitude high-speed cruise. According to the flight characteristics of different stages, the trajectory can be divided into different sections. Typical sections include:

(1) Program flight section

At the initial stage of takeoff and climb, the aircraft can fly according to a certain law of trajectory parameters. The flight program can construct different modes according to different parameters, such as the change of angle of attack α , the law of altitude H, etc. A typical variation law according to the trajectory inclination γ is: d γ d t = { C 1 0 < t < t 1 − C 2 t 1 ≤ t < t 2 . (3)

Among them, C ₁ and C ₂ can be taken as constants. Under this law, the aircraft takes off from the horizontal state, gradually decreases after reaching the maximum trajectory inclination, and finally turns into the level flight state.

If the height change rate is taken as the parameter, set the height change rate as a function of time, that is: d H d t = f ( t ) . (4) f ( t ) can be a constant value or the law of time. When the climbing ability is insufficient or you want to obtain a large acceleration rate, it can fly at constant altitude and the climb rate is zero, that is: d H d t = 0. (5)

(2) Variable acceleration flight

According to the performance of the engine, the acceleration rate V ˙ is taken as the control variable in the climbing process. The higher the thrust of the engine, the greater the acceleration rate that can be achieved, otherwise the acceleration rate is reduced. Acceleration rate as a function of time is achieved for a specific flight section: d V d t = f ( t , H , M a ) . (6)

If the acceleration rate is set to be constant, that is, a constant acceleration rate climb, that is: d V d t = C . (7)

(3) Isodynamic pressure flight

Isodynamic pressure flight is a common flight mode of aircraft, which can coordinate between acceleration rate and climb rate under the constraint of structural load. That is: { d Q d t = 0 Q 0 = C . (8)

(4) Cruise flight

When the aircraft reaches the predetermined cruise flight state, the flight altitude and speed remain constant. it is necessary to control the engine thrust through speed feedback and the balance relationship between aerodynamic force and moment to realize cruise flight. During cruise flight, the following requirements are met: { d H d t = 0 d V d t = 0 H = H C V = V C . (9)

In the above formula, H C is the cruise altitude and V C is the cruise speed.

(5) Transition Process Control

During flight, there will be differences in parameters between different sections. During the section conversion, the parameter PID feedback control is used to realize the smooth transition. For example, when transitioning from the constant acceleration phase to constant dynamic pressure flight, take Q ∗ as the expected dynamic pressure value by controlling the change of trajectory inclination Δ γ adjust the dynamic pressure. Construct the following model: Δ γ ∗ = K 1 ( Q − Q ∗ ) + K 2 Q ˙ + K 3 ∫ ( Q − Q ∗ ) d t . (10)

(6) Engine thrust

The main flight processes of horizontal take-off high-speed aircraft include ground take-off acceleration climb, constant speed cruise, return and other processes. In the climbing process, it is expected to climb at a large acceleration, and the engine works according to the maximum state, including: P = P M a x ( H , M a , Q , ε ) . (11) Where H, Ma, Q, and ε are flight altitude, Mach number, dynamic pressure, and fuel gas ratio, respectively, and P _Max is the thrust value under the maximum condition of the engine.

During cruise flight, the engine works in a throttling state, and the change in thrust can be calculated according to the feedback of a predetermined speed M a ∗ , and then the thrust can be adjusted through fuel supply. The thrust adjustment of the cruise section adopts the following form: Δ P = K 1 ( M a − M a ∗ ) + K 2 M a . (12)

(7) Constraints

During the flight, the flight profile, trajectory parameters, attitude angle, etc. will change. According to the design scheme, during the trajectory calculation, it is necessary to restrict the variation range of multiple parameters, mainly including:

Attack angle constraint: α ∈ [ α M i n , α M a x ] .

Dynamic pressure constraint: Q ∈ 1 2 ρ V 2 ∈ [ Q M i n , Q M a x ] .

Flight profile constraints: M a ≤ M a M a x ; H ≤ H M a x .

Overload restraint: N y ≤ N y M a x ; N z ≤ N z M a x .

Aerodynamic thermal restraint (Jia and Yan, 2015): q ˙ ≤ C 1 R d ( ρ ρ 0 ) 0.5 ( V V C ) 3.15 ≤ q ˙ M a x . (13)

The above constraints affect each other, so they are balanced according to certain strategies during flight.

For the flight process of horizontal take-off and landing on a high-speed cruise, the flight trajectories of different flight sections can be constructed. According to the model characteristics of different flight sections, the parameters affecting the flight process are extracted, and the parameter selection needs to be analyzed from the aspects of simplicity and sensitivity. Taking the longitudinal plane flight process as an example, this paper uses a typical four-section trajectory model for analysis. Since the range and flight time are mainly related to the climb and cruise process, the fuel threshold required for the return process is set in the calculation, and the return landing process is no longer compared. The main parameters are listed in Table 1, and the trajectory is shown in Figure 1.

The typical flight sections described above have a total of 11 parameters. By changing the parameter values and combining the constraints, different flight trajectories can be obtained. Obviously, in the flight profile, there are many parameters affecting the flight process, and there is a complex mutual coupling relationship between them.

Analyzing all parameters would greatly increase the difficulty of analyzing and optimizing the design. In practice, different parameters have different effects on flight trajectory. Through the analysis of a typical trajectory, four parameters are selected as the main variables for the trajectory analysis and optimization design for the flight mission with a certain cruise altitude and speed, including the flight time of the acceleration section T _A, the acceleration section acceleration V ˙ B S e t , the acceleration section end speed V _BEnd, and the dynamic pressure Q _Cset in the climb section.

The trajectory calculation model is a system of differential equations, and the fourth-order Runge-Kutta method is adopted for the differential equations. Let the initial value problem be expressed as follows: y ′ = f ( t , y ) , y ( t 0 ) = y 0 , y n + 1 = y n + h 6 ( k 1 + 2 k 2 + 2 k 3 + k 4 ) . (14)

Here: k 1 = f ( t n , y n ) , k 2 = f ( t n + h 2 , y n + h 2 k 1 ) , k 3 = f ( t n + h 2 , y n + h 2 k 2 ) , k 4 = f ( t n + h , y n + h k 3 ) .

The trajectory is solved by integrating on the time axis.

The BP neural network is a nonlinear parameter modeling method. Its most obvious feature lies in the error back-propagation learning algorithm it adopts, and it can adjust the weight coefficients of each layer network in the model in real time through continuous learning. When the total weight and fuel are constant, the variables of trajectory analysis and optimal design are used as input values, and the range R D and flight time T D are output values to establish a prediction model for the overall parameters of the trajectory. Since the input data of the neural network is given 4 parameters, the input layer has four nodes. The hidden layer is 1, the number of neurons is 8, and the output layer has two nodes. The adopted neural network structure is shown in Figure 2.

The input of the h neuron in the hidden layer is: a h = ∑ i = 1 4 ω i h x i , (15) where ω i h represents the weight of the i input neuron in the input layer to the h neuron in the hidden layer.

The activation function passing through the hidden layer is the tansig function, and the expression is: f ( x ) = 2 1 + e − 2 x − 1. (16)

Thus: b h = f ( α h − γ h ) , (17) where γ h represents the threshold of the h neuron in the hidden layer.

The input of the j neuron in the output layer is: β j = ∑ h = 1 8 v h j b h , (18) where v h j represents the weight from the h neuron in the hidden layer to the j output in the output layer.

The activation function of the output layer is the purelin function, and the expression is: f ( x ) = x . (19)

Thus, the output of the neural network is: y j ∼ = f ( β j − θ j ) , (20) where θ j represents the threshold of the j neuron in the output layer.

Establish loss function: J = 1 2 ∑ j = 1 2 ( y j ∼ − y j ) 2 , (21) where y j is the target output and y j ∼ is the output of the neural network.

By optimizing the input weights of neurons in each layer to minimize the loss function, the output of the neural network is close to the target output as much as possible, and the training model is obtained. Finally, the training model is used for prediction.

Based on the establishment of a parametric prediction model, optimization analysis can be carried out. There are many optimization design methods. Among them, the genetic algorithm, as a global optimization design method, has better optimization accuracy for nonlinear high-dimensional functions. In this paper, a genetic algorithm is used to optimize the neural network.

Through the combination of flight sections and parametric modeling, the trajectory optimization of the aircraft is transformed into the established neural network model and the process of optimization. The process of simulation calculation and modeling is as follows:

Step 1

Determine the model’s input and output parameters and sample points

According to the flight section analysis of the aircraft, taking the four main parameters that affect the flight section as the input and the range R _D and flight time T _D as the output, the functional relationship is established as follows: R D = f R ( T A , V ˙ B S e t , V B E n d , Q C S e t ) , T D = f T ( T A , V ˙ B S e t , V B E n d , Q C S e t ) . (22)

According to the working envelope of the aircraft, the analysis sample points for modeling are determined through experimental design or parameter combination. For the four parameter combinations in this paper, a total of 420 sample points are taken.

Step 2

Trajectory calculation of sample points

For the sample points, carry out the trajectory calculation in the flight process according to the trajectory calculation model established in Section 2, and obtain the sample values of the range R _D and flight time T _D . { R D 0 = f R 0 ( T A , V ˙ B S e t , V B E n d , Q C S e t ) R D 1 = f R 1 ( T A , V ˙ B S e t , V B E n d , Q C S e t ) ...... R D m = f R m ( T A , V ˙ B S e t , V B E n d , Q C S e t ) { T D 0 = f T 0 ( T A , V ˙ B S e t , V B E n d , Q C S e t ) T D 1 = f T 1 ( T A , V ˙ B S e t , V B E n d , Q C S e t ) ...... T D m = f T m ( T A , V ˙ B S e t , V B E n d , Q C S e t ) . (23)

Step 3

Establish neural network model based on sample points

Based on the trajectory calculation results of sample points, the neural network is trained to obtain the functional model between section parameters and range R _D and time T _D. On this basis, the established neural network model is used to predict the random sample points in the flight envelope. After comparing with the trajectory calculation results, the feasibility and accuracy of the model are analyzed.

Step 4

Model optimization

According to the established neural network model, the optimization of the neural network is carried out by using a genetic algorithm with flight range R _D as the optimization objective function. The trajectory calculation is carried out by using the segment parameter value corresponding to the best advantage obtained by optimization. The difference between the optimized value and the calculated value of range and time under the same segment parameters is compared, the feasibility of the optimization result is evaluated, and the characteristics of the optimized trajectory are analyzed.

In the calculation process, the neural network modeling and optimization algorithm parameters can be adjusted according to the verification of the model. The overall modeling and calculation process is shown in Figure 3.

For parameters T _A = 160s, speed V _BEnd = Ma1.8, and climb dynamic pressure Q _Cset = 40kPa, four different climb rates are used to calculate flight paths. The results are compared in Figure 4.

Before 160 s, the aircraft climbs according to the law of trajectory inclination. The flight Mach number continues to increase, and the dynamic pressure first increases and then decreases. The change is related to the inclination design in the program section. When compared with the trajectory of the climb section in Figure 4A, combined with the analysis of Mach number and dynamic pressure change, when the acceleration is V ˙ B S e t 0.5 m/s², the flight speed of Section 2 increases slowly. Under the maximum thrust of the engine, the climb rate of the aircraft is high, that is, the increase rate of height is large, so the dynamic pressure decreases rapidly in the initial stage, as shown in Figure 4C. Due to the low acceleration rate, when the flight time is 602.4 s, the speed reaches Ma1.8 and turns to Section 3 dynamic pressure flight. In this process, maintain a low dynamic pressure of about 20–25 kPa. Due to the long flight time of Section 2, the overall acceleration and climb time increases significantly, and the aircraft enters the cruise flight in 1300 s.

When the acceleration V ˙ B S e t of the constant acceleration section increases to 1 m/s, the climb rate of Section 2 decreases, the slope of the Mach number curve increases (Figure 4B), and the corresponding dynamic pressure also increases. Through calculation, the flight speed reaches Ma1.8 when the time is 377.9 s, and it turns to Section 3 constant dynamic pressure climb. Compared with the condition of acceleration of 0.5 m/s², the overall climbing time is significantly reduced, and it enters the cruise flight state at 1122 s.

When the acceleration of Section 2 is further increased, the climb rate of the aircraft is reduced under a certain thrust, and the kinetic energy is increased rapidly by reducing the increasing trend of potential energy. As shown in the trajectory curve in Figure 5A, when V ˙ B S e t = 2 m / s 2 , the constant acceleration section is approximately level flight, and when the acceleration is further increased to 4 m/s², a local dive is required to achieve a rapid increase in speed. From the change process of Mach number, when V ˙ B S e t = 2 m / s 2 , the flight speed reaches Ma1.8 at277.9 s and turns into a constant dynamic pressure flight section; when V ˙ B S e t = 4 m / s 2 , the flight speed is about Ma1.57 at 217.1 s, but the flight dynamic pressure has exceeded the set maximum dynamic pressure constraint value, so the flight speed is directly transfer to Section 3. From the perspective of dynamic pressure changes, after the acceleration rate exceeds 2 m/s², the dynamic pressure of Section 2 increases, and the maximum dynamic pressure exceeds the preset dynamic pressure value of Section 3. When transitioning to Section 3, the climb rate of the aircraft increases and the acceleration rate decreases until the dynamic pressure is restored. Decrease to the preset dynamic pressure value of Section 3, and then maintain the isodynamic pressure to fly. Because the acceleration time is shorter when the acceleration rate is large, the time to enter the cruise flight is also relatively early.

For parameters T _A = 160s, speed V _BEnd = Ma1.8, and acceleration V ˙ B S e t = 2 m / s 2 , four different climb rates are used to calculate flight paths. The results are compared in Figure 4. The flight trajectory under four groups of dynamic pressure Q C S e t is calculated.

From the change of parameters in Figure 5A, B, the flight dynamic pressure increases, the speed increases faster and the climbing time decreases. When Q C S e t = 25 kpa , the acceleration process is the longest. Due to the low dynamic pressure flight, the engine thrust is low and the climbing process is slow. It takes about 2693 s to reach the predetermined cruise flight state. With the increase of flight dynamic pressure, the engine thrust increases, and the climbing speed of the aircraft also increases. When the dynamic pressure is 60 kPa, the predetermined cruise parameters can be reached in 650 s, and the aircraft will turn to Section 4. From the trajectory in Figure 5A, in addition to the difference in climbing time and distance due to the cruise dynamic pressure setting of about 30 kPa and the climbing trajectory according to the dynamic pressure of 25 kPa, the maximum altitude has been higher than the cruise altitude before turning into cruise flight, and the flight altitude needs to be reduced when turning into Section 4. When the dynamic pressure during climbing is greater than 40 kPa, it needs to transition to the cruise altitude through further climbing because it is greater than the set cruise dynamic pressure. Under the given flight strategy, the dynamic pressure of Section 1 is first high and then low (Figure 5C). The dynamic pressure has achieved a good transition in the process of change, and there is no serious parameter overshoot and fluctuation. Corresponding to the changes in trajectory and Mach number, when turning to cruise flight, the dynamic pressure transits from the climb phase to the cruise dynamic pressure.

According to the ballistic simulation of the above typical state, the range and flight time under different combinations of parameters are extracted, which are compared with Figure 6.

From Figure 6A, the dynamic pressure Q C S e t of Section 3 has a great impact on the range. The range is significantly smaller under the low dynamic pressure of 25 kPa, which is due to the long climbing time and more total fuel consumption. When the dynamic pressure value of Section 3 is about 40 kPa, the range is the largest. If the dynamic pressure value is further increased, the range will be slightly reduced. From the influence of acceleration V ˙ B S e t , the influence under different dynamic pressures is different. When the flight dynamic pressure of Section 3 is 25 kPa, the range with an acceleration rate of 0.5 m/s² is the smallest, while the range with an acceleration rate of 2 m/s² and 4 m/s² is the largest. When the dynamic pressure is 40kPa, the range with an acceleration rate of 2 m/s² is the largest and the range with an acceleration rate of 1 m/s² is the smallest, which indicates that there is an interactive relationship between the design parameters of the flight section.

From Figure 6B, the dynamic pressure of Section 3 has a monotonic effect on the flight time. As the dynamic pressure Q C S e t increases, the flight time decreases. Compared with the influence of dynamic pressure, the influence of acceleration rate on flight time is relatively small. When the dynamic pressure is 60 kPa, the difference in flight time under different acceleration rates is less than 70 s.

According to the calculation results of sample points, a neural network prediction model corresponding to four parameters, flight range, and flight time is constructed. In the process of parameter construction, through the optimization of neural network model parameters, the main parameters are as follows:

1) Number of hidden layers is 1

In order to better test the prediction effect of the algorithm model, the experimental scheme of cross validation is adopted. The data is randomly divided into four subsets, which are recorded as S ₁, S ₂, S _3, and S ₄. Select S ₁, S ₂, S _3, and S ₄ as the test sets, respectively, and select the other three subsets as the training set to establish the model. The experimental and verification scheme are shown in Table 2.

The neural network is trained according to the test and verification scheme given in Table 2. The R _D and T _D training and prediction results of the two sets of trials are compared in Figure 7 and Figure 8. From the results, the training model’s error in the sample prediction value of the flight range is mainly concentrated in ±10 km, and the value error of the flight time prediction value is mainly between ±5 s. The error value of individual test sample points is relatively large, but the relative error is not high, indicating that the neural network has better accuracy.

The random state test is carried out for the established prediction model. Table 3 shows the range and flight time calculated by using neural network and trajectory model, respectively, under the randomly selected 10 groups of parameter combinations. From the comparison of results, the maximum error between the range prediction value R _D(BP) and the trajectory calculation value R _D(Tra) is in condition 4. The relative error of the two algorithms is less than 0.82%, and the error of other conditions is less than 0.3%. The maximum relative error between the predicted value T _D(BP) of flight time and the calculated value T _D(Tra) of trajectory is within 0.45%. It can be seen that the neural network prediction model has high accuracy.

The aircraft’s range is an important indicator of the overall design. In the trajectory design, overload, attack angle, dynamic pressure, etc. have been reflected in the flight model as constraints, so the range R D is used as the objective function of trajectory optimization.

A genetic algorithm is used to optimize the overall parameters of the established neural network prediction model. The main parameters of the algorithm refer to the values in the study by Cheng and Wang (2011), and the settings are listed in Table 4:

In order to test the influence of genetic algorithm parameters, the group sizes of 50, 100, 150, 200, and 250 are taken, and the optimal value of R _D obtained by optimization varies from 4991.40 to 4992.26 km; Take five groups of mutation probability of 0.05, 0.10, 0.15, 0.20, and 0.25. The variation range of R _D is 4991.40–4992.78 km, and its relative variation value is very small. For the training model, the change of algorithm parameters is not sensitive to the optimization results, so it can be carried out according to the parameter values in Table 4.

Based on the parameter settings in Table 4, a total of 138 steps are iterated, and the calculated results are listed in Table 5. From the optimization results, the flight pressure is close to 40 kPa, which is similar to the ballistic characteristics analysis results in Section 4.

The optimal value in Table 5 is used as the input parameter for trajectory calculation. Under this condition, the flight time of the aircraft is 3247.60 s and the range is 4981.15 km. The result of the trajectory calculation is basically consistent with the time prediction value in Table 4. The range value is slightly smaller, and the relative error is 0.2058%. Comparing the overall parameters of the sample points, the optimized results, and the optimal point of trajectory calculation in Figure 9, it can be seen that the optimized result range R D is the best, and the flight time T D is at the middle level of the sample points.

The trajectory and parameter changes in the optimized state are shown in Figure 10. From the analysis of the flight process, since the end speed of Section 2 is Ma1.205, the proportion of this section in the climb process is relatively small, and it is transferred to Section 3 isodynamic flight when the flight altitude is about 8.7 km. During the whole climbing process, the flight speed continued to increase, and the acceleration in the isodynamic pressure section was large. After the change of dynamic pressure, the maximum value is about 55 kPa, which does not exceed the upper and lower limits of constraints, and the parameter changes are within a reasonable range.