Nomenclature
-
${h_0}$
-
atmospheric density elevation [m]
-
${R_e}$
-
radius of the earth [m]
-
$\dot Q$
-
heat flow [W]
- q
-
dynamic pressure [Pa]
-
${N_y}$
-
overload [g]
-
${S_a}$
-
reference area [
${{\rm{m}}^2}$
] - L
-
lift force [N]
- D
-
drag force [N]
-
${C_D}$
-
lift coefficient
-
${C_L}$
-
drag coefficient
- Ma
-
Mach number
- RMSE
-
root mean square error
- MAE
-
mean absolute error
- MAPE
-
mean absolute percentage error
- ABC
-
Artificial Bee Colony
- MLP
-
multilayer Perceptrona
Greek symbols
-
$\alpha $
-
angle-of-attack [
$ \circ $
] -
$\sigma $
-
angle-of-bank [
$ \circ $
] -
$\gamma $
-
flight path angle [
$ \circ $
] -
$\psi $
-
heading angle [
$ \circ $
] -
$\theta $
-
longitude [
$ \circ $
] -
$\phi $
-
latitude [
$ \circ $
] -
$\rho $
-
atmospheric density at sea level [
${\rm{kg/}}{{\rm{m}}^3}$
] -
$\mu $
-
random number ranging from 0 to 1
1.0 Introduction
A hypersonic glide vehicle is a lifting flight that can travel in near space with high speed, long range, and good penetration ability [Reference Liang, Ren, Li and Chen1, Reference Yu, Dong, Li, Ren and Lv2]. When a hypersonic glide vehicle is used as an offensive weapon, its irregular orbital flight significantly increases the difficulty of ballistic prediction and interception by existing anti-missile air defense systems [Reference Zhao and Zhou3].
The trajectory optimisation of a hypersonic gliding vehicle is a nonlinear optimal control problem with multiple complex constraints. In the recent years, many researchers have conducted extensive research in this area; furthermore, this problem can solved using two types of methods: direct and indirect methods. The indirect method is more accurate; however, the derivation process is complex, and the algorithm is sensitive to the initial values and poorly convergent. The direct method transforms the optimal control problem into a nonlinear programming problem by directly discretising control and/or state quantities. Pseudospectral method, a direct method, is widely adopted in the aerospace field. Zhao et al. [Reference Zhao, Huang, Li and Yang4] proposed a set of re-entrant trajectory planning methods under complex constraints using the Gaussian pseudospectral method. The method was described in [5] for computational optimal guidance and control using adaptive Gaussian quadrature collocation and sparse nonlinear programming.
However, as the anti-missile weapon technology continues to improve and develop with time, the survivability of individual craft decreases. Consequently, hypersonic vehicles can be employed in a saturation attack style of operation, where the combined effectiveness of the cluster is effectively increased by simultaneous strikes from multiple angles or batches. Research on multibomb synergy algorithms has been underway since 2005. Tactical missiles have been studied to design time-constrained guidance laws based on the small heading angle assumption and optimal control theory to achieve multibomb time-of-flight synergy [Reference Gregg6, Reference Kumar and Ghose7]. Furthermore, an optimal single-bullet guidance law with a fixed total time constraint was designed to achieve the simultaneous arrival of multiple bombs at the target point [Reference Jeon, Lee and Tahk8].
In contrast to multi-bomb and multi-unmanned aerial vehicle (UAV) coordination, hypersonic glide vehicles are under-actuated systems with limited control and a more complex flight environment, making it impossible to directly apply multi-bomb coordination methods. Hypersonic vehicles travel at speeds greater than Mach 5 at a height range of 20–100km above the ground; therefore, their trajectory differs from that of conventional aircraft. Generally, the trajectories of hypersonic vehicles can be divided into two types: cruising and gliding. Both types of flights require a booster section to bring the vehicle to a predetermined altitude, followed by powered cruise flight or unpowered glide flight. This study focused on the trajectory of hypersonic gliding vehicles. Hypersonic gliding vehicles have a wide range of altitude changes, unpredictable flight paths, fast flight speeds and relatively long ranges. Moreover, they have significant advantages in terms of long-range rapid precision strikes and power delivery.
To achieve this operational style, many scholars have conducted research on collaborative guidance methods for hypersonic vehicles, primarily through a method known as time-constrained guidance. Zhao et al. proposed a cooperative guidance method based on the premise of normal communication between [Reference Zhao and Zhou9, Reference Zhao, Zhou and Dong10], where a cooperative attack on the target is achieved by adopting different coordination variables and control methods. Liu et al. studied the cooperative guidance problem of multiple inferior missiles intercepting a hypersonic target with a specific impact angle constraint in a two-dimensional plane, and designed the fractional-order sliding mode guidance law for intercepting hypersonic vehicles [Reference Liu, Yan, Zhang, Liu and Yan11, Reference Liu, Yan, Dai, Yan and Xin12]. Han et al. proposed a three-dimensional (3-D) guidance law that considers the terminal angle constraint against manoeuvering targets. Each vehicle achieves a simultaneous hit on the target using a time-constrained guidance law based on a predetermined desired flight time [Reference Han, Hu, Shin, Tsourdos and Xin13, Reference Han, Hu and Xin14]. An and Xu et al. [Reference An, Guo, Huang and Xu15–Reference An, Guo, Huang and Xu17] studied a distributed time cooperative guidance method based on multi-agent cooperative control theory and obtained a robust tracking control scheme for the hypersonic gliding leap trajectory.
Although the aforementioned guidance methods can achieve good temporal synergy, they ignore the aerodynamic characteristics and are only suitable for solving the problem of the end guidance segment [Reference Dalle, Torrez, Driscoll, Bolender and Bowcutt18, Reference Franco, Rivas and Valenzuela19]. However, the end-guidance segment of a hypersonic glide vehicle only occupies a small fraction of the flight time; therefore, it is not sufficient to carry out collaborative planning only at the end of the flight [Reference Ng, Sridha and Grabbe20]. Therefore, it is of great importance to study collaborative trajectory optimisation methods applicable to hypersonic gliding vehicles.
To solve the problem of cooperative trajectory optimisation of hypersonic gliding vehicles, this study first obtained a reasonable cooperative time as the constraint of the optimisation problem through multilayer Perceptrona (MLP) neural network training. Additionally, the neural network determines the optimal initial heading angle of each aircraft to reduce unnecessary manoeuvering performance consumption in the flight mission. Furthermore, to obtain satisfactory results, the attack angle and bank angle commands are set as functions of speed, and some parameters of the functions are optimised through Artificial Bee Colony(ABC) algorithm. The key contributions of this study are summarised as follows:
-
(1) A collaborative trajectory planning model for hypersonic gliding vehicles was developed including dynamical differential equations and associated constraints.
-
(2) To obtain the effective collaborative time constraints, an MLP model was trained and adopted in a time decision module, whose input is a four-dimensional vector selected according to the trajectory characteristics. The MLP also determines the optimal initial heading angle of each aircraft to reduce unnecessary manoeuvering performance consumption in the flight mission.
-
(3) A three-dimensional trajectory optimisation strategy is designed. The angles of attack and banks are designed as velocity-dependent profiles. The adjustable parameters are optimised using the ABC algorithm to improve the accuracy of the coordinated trajectory optimisation under fixed time constraints.
The remainder of this paper is organised as follows. Section 2 presents the modeling process for the optimisation of cooperative trajectories of hypersonic glide vehicles. Section 3 describes the method to determine the cooperative flight time and optimal initial heading angle by MLP and the setting and optimisation method of the angle-of-attack and bank. Section 4 presents the results of the optimisation and comparative analyses. Section 5 presents concluding remarks.
2.0 Problem formulation
2.1 Equations of motion
Assuming that Earth is a rotating sphere, the unpowered three-dimensional motion equations of each hypersonic gliding vehicle flying in concert can be expressed by Equation (1).
In Equation (1), i indicates the serial number of the vehicle in the cluster, r is the geocentric distance,
$\theta $
and
$\phi $
are the longitude and latitude, respectively, corresponding to the current position of the aircraft. The flight-path and velocity-heading angles are represented by
$\gamma $
and
$\psi $
, respectively.
$\sigma $
is the bank angle and
${\omega _e}$
is the rotation velocity of the earth, which is assumed to be a fixed value.
$g = \mu /{r^2}$
is the acceleration due to gravity and
$\mu $
is the Earth’s gravitational constant.
\begin{align} \left\{ {\begin{array}{l}{\dfrac{{d{r_i}}}{{dt}} = } {{v_i}\sin {\gamma _i}}\\[12pt] {\dfrac{{d{\theta _i}}}{{dt}} = } {\dfrac{{{v_i}\cos {\gamma _i}\sin {\psi _i}}}{{{r_i}\cos {\phi _i}}}}\\[12pt] {\dfrac{{d{\phi _i}}}{{dt}} = } {\dfrac{{{v_i}\cos {\gamma _i}\cos {\psi _i}}}{{{r_i}}}}\\[12pt] {\dfrac{{d{v_i}}}{{dt}} = } { - {D_i} - {g_i}\sin {\gamma _i} + {\omega _e}^2r\cos {\phi _i}(\sin {\gamma _i}\cos {\phi _i} - \cos {\gamma _i}\sin {\phi _i}\cos {\psi _i})}\\[12pt] {\dfrac{{d{\gamma _i}}}{{dt}} = } {\dfrac{{{L_i}\cos {\sigma _i}}}{{{v_i}}} - \dfrac{{{g_i}}}{{{v_i}}}\cos {\gamma _i} + \dfrac{{{v_i}}}{{{r_i}}}\cos {\gamma _i} + 2{\omega _e}\cos {\phi _i}\sin {\psi _i}}\\[12pt] {} \quad { + \dfrac{{{\omega _e}^2{r_i}}}{{{v_i}}}\cos {\phi _i}(\cos {\gamma _i}\cos {\phi _i} + \sin {\gamma _i}\cos {\psi _i}\sin {\phi _i})}\\[12pt] {\dfrac{{d{\psi _i}}}{{dt}} = } {\dfrac{{{L_i}\sin {\sigma _i}}}{{{v_i}\cos {\gamma _i}}} + \dfrac{{{v_i}}}{{{r_i}}}\cos {\gamma _i}\sin {\psi _i}\tan {\phi _i} - 2{\omega _e}(\tan {\gamma _i}\cos {\psi _i}\cos {\phi _i} - \sin {\phi _i})}\\[12pt] \quad {} { + \dfrac{{{\omega _e}^2}}{{{v_i}\cos {\gamma _i}}}\sin {\psi _i}\sin {\phi _i}\cos {\phi _i}}\end{array}} \right.\end{align}
L and D are the lift acceleration and drag acceleration, respectively, which can be calculated using the following formula:
\begin{align} \left\{ {\begin{array}{l}{D = \dfrac{1}{{2m}}\rho {V^2}{S_a}{C_D}(\alpha ,Ma)}\\[12pt] {L = \dfrac{1}{{2m}}\rho {V^2}{S_a}{C_L}(\alpha ,Ma)}\end{array}} \right.\end{align}
where
${S_a}$
is the dimensional reference area of the vehicle and
${C_L}$
and
${C_D}$
are the lift and drag coefficients, respectively, which are functions of the attack angle
$\alpha $
and Mach number Ma.
$\rho $
is the atmospheric density, which can be calculated as:
\begin{align} \rho = {\rho _0}{e^{((r - {R_e})/{h_0})}}\end{align}
where
${\rho _0}$
is the atmospheric density at sea level (
$1.752kg/{m^3}$
) and
${h_0}$
is the elevation at which the density is being evaluated (6700m).
2.2 Constraints
2.2.1 Control variant constraints
When optimising the glide phase trajectory, the two control variables of each vehicle are the attack angle
${\alpha _i}$
and bank angle
${\sigma _i}$
. In this stage, both the attack and bank angles should have an appropriate upper and lower boundary to ensure that the vehicle has sufficient distance to fly both longitudinally and laterally. For example, the vehicle will not have sufficient lift to achieve the desired range if the attack angle is excessively small. Furthermore, an excessively large attack angle will increase both the lift and drag forces, which will result in excessively violent ballistic jumps and loss of vehicle speed. The control variables are set as velocity-dependent variation curves. Their amplitude limits can be expressed by Equation (4) and the specific design and optimisation strategy of the control volume are described in detail in the following section.
\begin{align} \left\{ {\begin{array}{c}{{\sigma _{\min }} \le \left| {{\sigma _i}} \right| \le {\sigma _{{\rm{max}}}}}\\[5pt] {{\alpha _{\min }} \le \left| {{\alpha _i}} \right| \le {\alpha _{{\rm{max}}}}}\end{array}} \right.\end{align}
2.2.2 Path constraints
Hypersonic gliding vehicles have high speed and large manoeuvering range. Their flight environment is complex, and the external factors and limitations of its own structure and material strength can affect the flight process [Reference Zhang, Yan, Huang, Che and Wang21]. Negative effects can be avoided throughout the flight of vehicles by adding path constraints. Typical path constraints include the heating rate
${\dot Q_i}$
, aerodynamic overload
${N_{yi}}$
and the dynamic pressure limit
${q_i}$
. Thus, the path constraints for the ith vehicle can be expressed by Equation (5).
\begin{align} \left\{ {\begin{array}{l}{{{\dot Q}_i} = {K_Q}\sqrt {{\rho _i}} {V_i}^{3.15} \le {{\dot Q}_{\max }}}\\[8pt] {{q_i} = \dfrac{1}{2}\rho {V_i}^2 \le {q_{\max }}}\\[12pt] {{N_{yi}} = \dfrac{{({C_L}\cos {\alpha _i} + {C_D}\sin {\alpha _i}){\rho _i}{V_i}^2{S_a}}}{{2m{g_0}}} \le {N_{\max }}}\end{array}} \right. \end{align}
where
${K_Q}$
is a constant with respect to the structural property of the gliding vehicle,
${V_i}$
and
$kg/{m^3}$
are the units of m/s and atmospheric density,
${\rho _i}$
, respectively.
2.2.3 Initial and terminal constraints
To accomplish the coordinated mission, all the vehicles must be able to satisfy the latitude, longitude and altitude constraints at the endpoint [Reference An, Guo, Xu and Huang22]. Additionally, the initial state of vehicles should be constrained. The initial and terminal constraints can be expressed as
\begin{align} \left\{ {\begin{array}{c}{{h_i}({t_0}) = {h_{0,i}}}\\[5pt] {{\theta _i}({t_0}) = {\theta _{0,i}}}\\[5pt] {{\phi _i}({t_0}) = {\phi _{0,i}}}\end{array}} \right. \end{align}
\begin{align} \left\{ {\begin{array}{c}{{h_i}({t_f}) = {h_{f,i}}}\\[5pt] {{\theta _i}({t_f}) = {\theta _{f,i}}}\\[5pt] {{\phi _i}({t_f}) = {\phi _{f,i}}}\end{array}} \right. \end{align}
where
${t_0}$
and
${t_f}$
denote the initial flight time and arrival time at the target point. For the mission requirement of k vehicle time synergy, the condition of Equation (8) needs to be satisfied:
\begin{align} {t_{f,1}} = {t_{f,2}} = \cdots = {t_{f,k}} = {t^*} \end{align}
where f represents the terminal condition,
$1,2 \cdots k$
represents the serial number of the vehicles, and
${t^*}$
represents the cooperative flight time.
3.0 Time coordinated trajectory optimisation strategy
This study proposes a trajectory planning strategy with a cooperative arrival time constraint to realise the mission requirements of multiple glide hypersonic vehicles. The overall optimisation process is illustrated in Fig. 1. First, the MLP model is trained and adopted to obtain the cooperative flight time and optimal initial heading angle of each aircraft. Then, both the angle-of-attack and bank are set in speed-dependent forms, and the adjustable parameters are optimised by the ABC algorithm to improve the cooperative flight accuracy of both cross-range and downrange.
Figure 1. Collaborative trajectory optimisation scheme.
3.1 Neural network model construction
Cooperative flight time is an extremely important metric when performing co-trajectory optimisation of hypersonic glide vehicles. Meanwhile, if the optimal initial heading angle for the vehicle can be selected during trajectory optimisation, the vehicle can reach the desired destination with less adjustments, which will reduce the complexity and difficulty of vehicle control. Therefore, the MLP model was trained to obtain the cooperative flight time and optimal initial heading angle simultaneously. The MLP structure is illustrated in Fig. 2.
Figure 2. Neural network structure diagram.
MLP neural network is one of the most widely used and rapidly developing artificial neural networks, and both theoretical research and practical applications have reached a high level. The sample quality can significantly affect the generalisation ability of neural networks; therefore, improving the training sample quality is an important method for improving the generalisation ability of neural networks. The inputs to the neural network should be selected based on ballistic characteristics and have a suitable relationship with the outputs.
In the cooperative trajectory optimisation task of multiple hypersonic glide vehicles, the initial and target positions of each vehicle are known, and the starting and ending positions of the glide vehicles have a mapping relationship with flight direction and time. However, the trajectories of hypersonic gliders are extremely irregular in the longitudinal direction; therefore, the sample quality will decline if altitude data are used as the MLP input. As a result, only the longitude and latitude data were used as the input of the neural network; that is, the downrange data of the aircraft were not directly used as the input of the neural network.
The samples for the neural network training were obtained through simulation calculations. First, the initial longitude, latitude, heading angle and flight time data of the vehicles were obtained by sampling. Second, the obtained data were substituted into the differential equations in the dynamics to solve for the end-position data of the vehicle. Subsequently, the input and output samples of the neural network were obtained by classifying all data as required. In this study, 53,760 sets of samples were generated for the neural network training. Figure 3 shows a part of the calculated sample data.
Figure 3. Sample acquisition process.
3.2 3-D trajectory optimisation based on ABC algorithm
3.2.1 Command design for angle-of-attack and bank
After determining the cooperative flight time using the MPL model, the cooperative arrival accuracy of the vehicles needs to be improved by optimising the angle-of-attack and bank. This subsection describes the optimisation methods for angle-of-attack and bank. In previous studies, the command curves of the angle-of-attack and bank were set as piecewise functions that varied with speed, and the speed values at the turning points were fixed. However, this type of command form cannot be adjusted adequately and cannot obtain the best results for all the tasks. Additionally, in this study, the angle-of-attack and bank reference profiles were set as functions of velocity; the difference is that some adjustable parameters were added. The optimal values are obtained using the ABC algorithm to satisfy the initial and terminal constraints by adjusting the angle-of-attack and bank commands. The angle-of-attack was set as Equation (9), taking the form of velocity-dependent segmentation functions.
\begin{align} \alpha = {\rm{10}} - {k_1}/({V_f}/{\rm{100}} - {V_{ref1}}/{\rm{100}}{)^2}*((V - {V_{ref2}})/{k_2}{)^2}\end{align}
In Equation (9),
${k_1}$
,
${k_2}$
,
${V_{ref1}}$
, and
${V_{ref2}}$
are the adjustable parameters. Similar to the attack angle curve, the setting of the bank angle command is related to the speed of the vehicles. The amplitude of the bank angle affects the longitudinal motion capability of the aircraft, while changes in the positive and negative signs affect the lateral motion capability. After determining the optimal initial heading angle and angle-of-attack command, the aircraft was able to reach the target point attachment in both crossrange and downrange. To further improve the precision of the coordinated arrival, the bank angle command is set in the form of a piecewise constant value, as shown in Equation (10).
\begin{align} \sigma = \left\{ {\begin{array}{c}{\begin{array}{l}{{\rm{0}}\quad {\rm{V}} \gt {{\rm{V}}_{\sigma 1}}}\\[5pt] {{\sigma _{ref1}}\quad {{\rm{V}}_{\sigma 2}} \le {\rm{V}} \le {{\rm{V}}_{\sigma 1}}}\\[5pt] {{\sigma _{ref2}}\quad {\rm{V}} \lt {{\rm{V}}_{\sigma 2}}}\end{array}}\end{array}} \right.\end{align}
${\sigma _{ref1}}$
,
${\sigma _{ref2}}$
,
${{\rm{V}}_{\sigma 1}}$
,
${{\rm{V}}_{\sigma 2}}$
are adjustable parameters and the optimal values are selected by the ABC algorithm.
3.2.2 ABC algorithm
The ABC algorithm is a swarm optimisation algorithm having the advantages of few control parameters, strong robustness and easy implementation [Reference Duan and Li23, Reference Wu, Wang, Li, Yao, Huang, Su and Yu24]. Furthermore, it has produced good results in many application fields. There are four sequentially realised phases in the algorithm: initialisation, employed bees, onlooker bees, and scout bee phases. In the initialisation phase, the dimensions of the optimisation problem m, number of nectar sources SN, and position vector of each nectar source
${{\bf{X}}_i} = [{x_{i1}},{x_{i2}}, \cdots ,{x_{im}}](i = 1,2, \cdots ,SN)$
are defined. The initial nectar source position was randomly generated according to Equation (11).
\begin{align} {x_{id}} = l{b_d} + (u{b_d} - l{b_d}) \cdot rand(0,1)\end{align}
where,
$i = 1,2, \cdots ,SN$
and
$d = 1,2, \cdots ,m$
.
${x_{id}}$
is the d-dimensional vector of
${{\bf{X}}_i}$
,
${\bf{ub}} = [u{b_1},u{b_2}, \cdots u{b_m}]$
and
${\bf{lb}} = [l{b_1},l{b_2}, \cdots u{b_m}]$
represent the upper and lower limits of the dimensional search space, respectively.
$rand(0,1)$
denotes a random number in [0,1].
To obtain a new nectar source
${v_i}$
, the employed bee at nectar source
${x_i}$
randomly selects another bee from the honey source
${x_k}$
to search the neighborhood and update the location. The search process for the employed bees can be expressed by Equation (12).
\begin{align} {v_{i,d}} = {x_{i,d}} + ({x_{i,d}} - {x_{k,d}}) \cdot rand(0,1)\end{align}
where
${v_{i,d}}$
is the d-dimensional vector of
${{\bf{v}}_i}$
and
$rand(0,1)$
is a random number in [0,1].
$k \in \left\{ {1,2, \cdots ,SN} \right\}$
and
$k \ne i$
is satisfied. Updating the new nectar source according to greedy selection and trail is the value of the discard counter, indicating that the quality of the honey source did not improve after multiple searches.
After a round of updating, the honey source information is shared with the onlooker bee, which selects the onlooker bee makes probability according to the quality of the nectar source. The probability of the second nectar source being selected by the onlooker bee is calculated using Equation (13).
\begin{align}{p_i} = \frac{{fi{t_i}}}{{\sum\limits_{j = 1}^{SN} fi{t_j}}}\end{align}
The fitness value
$fi{t_j}$
was calculated using Equation (14):
\begin{align}\left\{ {\begin{array}{c}{\dfrac{1}{{1 + {f_i}}},{f_i} \geqslant 0}\\[12pt] {1 + \left| {{f_i}} \right|,{f_i} \lt 0}\end{array}} \right.\end{align}
where
${f_i}$
is the evaluation value of the ith nectar source calculated from the objective function of the problem to be solved. The working diagram of the basic ABC algorithm is shown in Fig. 4.
Figure 4. Optimisation flow chart of ABC algorithm.
3.2.3 ABC algorithm based optimisation strategy
The collaborative trajectory optimisation process based on ABC algorithm under fixed time constraints is shown in Fig. 5, in which the subsystem is trajectory calculation module. After initialising the algorithm, the subsystem calculates the initial nectar position to form the current optimal solution. Then, the nectar source position in the stage of employed and onlooker bees determines whether there are scout bees. If there is a scout bee, the optimisation is transferred to the scout bee stage. Otherwise, it is judged whether the optimisation end conditions are met, and if the optimisation results do not meet the index requirements, the next iteration is required.
Figure 5. 3-D trajectory optimisation flow chart.
4.0 Results and discussion
The same type of hypersonic glide vehicle with the same structural and aerodynamic parameters is assumed to be involved in the collaborative planning. In this study, the general and aerodynamic parameters of CAV-H [Reference Phillips25], a general aviation vehicle solution, was used for the simulation and validation. The vehicle is rocket-boosted to near space and then begins its long-range, unpowered glide through the atmosphere. The mass of the vehicle was 907kg, the reference area was
$0.484\,{m^2}$
and the maximum lift-drag ratio was 3.5. Unless specified otherwise, all the numerical results presented in this paper were generated using a computer with an Intel i7-11800H processor (2.3 GHz, 8 cores, 16 threads), 16.0 GB memory, and Windows operating system.
4.1 Analysis of collaborative time-of-flight and optimal initial heading angle prediction models
In this study, we used an artificial neural network with four hidden layers, and the training function of all layers was a Rectified Linear Unit. The number of nodes in the four hidden layers was 32, 64, 128 and 64. After sampling 53,760 groups of data according to the method in Section 3.1, 44,800 groups (83.33%) were randomly selected as the training database and 8,960 groups (16.67%) were selected as the validation database. The results of the predictive model training are shown in Figs. 6–9.
Figure 6. Training loss changes during training.
Figure 7. Root mean square changes during training.
Figure 8. Mean absolute error changes during training.
Figure 9. Mean absolute percentage error changes during training.
In the figures, some variables are used to describe the training process and performance of the neural network. The loss function measures the degree of inconsistency between the true and predicted values, RMSE, MAE and MAPE are denoted as the root mean square error, mean absolute error, and mean absolute percentage error, respectively, which can be calculated using Equations (15) to (17).
\begin{align} RMSE = \sqrt {\frac{1}{N}\sum\limits_{i = 1}^N {{({y_i} - {{\tilde y}_i})}^2}} \end{align}
\begin{align} MAE = \frac{1}{N}\sum\limits_{i = 1}^N \left| {{y_i} - {{\tilde y}_i}} \right|\end{align}
\begin{align} MAPE = \frac{1}{N}\sum\limits_{i = 1}^N \left| {\frac{{{y_i} - {{\tilde y}_i}}}{{{y_i}}}} \right|\end{align}
The training process for the neural network is illustrated in Figs. 6 and 7. In Fig. 6, the loss function values are close to zero at the final stage of training, indicating that good convergence was obtained. It can be observed in Fig. 7 that when the iteration number of the validation group is 399, the root mean square error is less than 0.05, thus meeting the application requirements in this paper.
The distributions of MAE and MAPE in the validation data when the neural network reached the required accuracy are shown in Figs. 8 and 9, indicating that the neural network proxy model has good fitting accuracy.
4.2 Collaborative trajectory optimisation simulation
Three hypersonic glide vehicles were used in the trajectory optimisation simulation experiment, and their initial conditions are listed in Table 1, respectively. The terminal constraints of the three vehicles were consistent, and the specific information is presented in Table 1.
Table 1. Initial state of each vehicle and target point state
Figure 10. Optimised collaborative trajectories.
Figure 11. Angle-of-attack comparison.
Figure 12. Angle-of-bank comparison.
Figure 13. Flight range comparison.
Figure 14. Final speed comparison.
The same collaborative trajectory optimisation problem was solved using the GPOPS toolbox based on the Gaussian pseudospectral method and the strategy proposed in this study to facilitate a comparison of their results. The optimisation strategy proposed in this study was denoted as M-ABC. In the ABC algorithm part of the optimisation strategy, the initial population size was set to 15, and the maximum number of iterations was set to 100. The coordinated trajectories optimised using the GPOPS toolbox and M-ABC are shown in Fig. 10.
It can be observed in Fig. 10 that although each vehicle can reach the target point accurately, the flight path optimised by the two strategies for each vehicle is different. As the optimisation strategy proposed in this study optimises the initial heading angle of the vehicle, the flight is oriented in a more optimal direction from near the starting point, diverging from the trajectory formed by the GPOPS toolbox. Simultaneously, different attack angles and bank angle commands give vehicles different performances throughout their trajectories. The strategy proposed in this paper sets the attack angle and the bank angle commands as speed-dependent functions, optimising only specific variables in the function, which will significantly reduce the complexity of the command curve resulting from the optimisation. A simpler command will result in a more easily executed trajectory, which is flatter and does not unnecessarily deplete the manoeuverability of the vehicle.
The angle-of-attack and bank command curves obtained using the two strategies are shown in Figs. 11 and 12. The angle-of-attack command mainly affects the longitudinal motion capability of the aircraft, whereas the amplitude of the bank angle affects the longitudinal motion capability of the aircraft, and changes in the positive and negative signs affect the lateral motion capability. The angle-of-attack and bank command curves optimised by GPOPS toolbox have poor smoothness, which is unfavourable for the design of aircraft control systems. Evidently, the attack angle and bank angle command curves obtained by M-ABC are better than the optimisation results of GPOPS toolbox in terms of smoothness. It will be easier to design an attitude-control algorithm for vehicles based on such instructions.
Figure 13 shows the comparison of the flight ranges of the vehicles obtained by the two optimisation strategies. Because the optimisation of the initial heading angle of the aircraft is considered in the optimisation scheme designed in this study, each vehicle can reach the target point within a relatively short range. The ranges of the three aircraft were shortened by 12.04%, 2.95% and 24.25%, with an average of 13.08%.
The manoeuverability of a vehicle is primarily reflected by its speed. Figure 14 shows the terminal speeds of the vehicles using the two optimisation strategies. It can be seen from the figure that the strategy proposed in this study can increase the terminal speeds of the three vehicles. Compared to the results of GPOPS toolbox, the final speeds of the three vehicles are increased by 423.3, 44.1 and 479.4m/s, respectively, with an average value of 315.6m/s.
5.0 Conclusions
This study proposed a cooperative trajectory optimisation method for hypersonic gliding vehicles with fixed time constraints. When establishing the problem model, the associated constraints were considered based on the differential equations of the vehicle dynamics. Moreover, an MLP artificial network was trained to obtain a reasonable cooperative time. Additionally, to avoid unnecessary range and speed consumption simultaneously, the MLP neural network can obtain the optimal initial heading angle for each vehicle. The 3-D trajectory optimisation problem under cooperative time is solved by optimising the speed-dependent angle-of-attack and bank based on the ABC algorithm. Numerical simulation with the CAV-H model showed that, compared to the GPOPS toolbox, the optimisation scheme proposed in this paper can obtain more smooth control commands, and each aircraft saves an average of 13.08% of the flight range, and the terminal speed is increased by 315.6 m/s, which can effectively improve the safety during flight and terminal manoeuver capability.
Acknowledgements
The authors appreciate the financial support from the Natural Science Basic Research Plan in Shaanxi Province of China (Grant No. 2020JC-19), National Natural Science Foundation of China (NSFC) (Grant No. 62173274), and Natural Science Foundation of Shaanxi Province (Grant No. 2020JQ-219).
