2.1 Overview of ChOA

The ChOA is a swarm-based optimisation algorithm inspired by the intelligence and sexual motivation of chimps in their crowd hunting behaviours. A chimp swarm includes four types of chimps: driver, barrier, chaser and attacker, with different skills necessary for a successful hunt. The attacker is the leader of the swarm. The other three types of chimps assist in hunting, and their social status decreases in the order of barrier, chaser and driver. The hunting process can be divided into diversification (exploration) and intensification (exploitation) stages. The exploration stage includes driving, blocking, and chasing, and the exploitation process is focused on attacking the prey. The prey is hunted in both stages. The processes of hunting and chasing the prey can be modelled using Eqs. (1) and (2).

(1) \begin{eqnarray}\boldsymbol{d}=\left| \boldsymbol{c}\cdot{\boldsymbol{x}_{prey}}(t)-\boldsymbol{m}\cdot{\boldsymbol{x}_{chimp}}(t) \right|\!,\end{eqnarray}

(2) \begin{eqnarray}{\boldsymbol{x}_{chimp}}(t+1)={\boldsymbol{x}_{prey}}(t)-\boldsymbol{a}\cdot\boldsymbol{d},\end{eqnarray}

where t is the number of the current iteration; $\boldsymbol{a}$ , $\boldsymbol{m}$ , and $\boldsymbol{c}$ are the coefficient vectors; $\boldsymbol{x}_{prey}$ represents the prey position; and $\boldsymbol{x}_{chimp}$ represents the chimp position vector. $\boldsymbol{a}$ , $\boldsymbol{m}$ , and $\boldsymbol{c}$ can be calculated as

(3) \begin{equation}\boldsymbol{a}=2\cdot{f}\cdot{\boldsymbol{r}_{1}}-{f,} \end{equation}

(4) \begin{equation}\boldsymbol{c}=2\cdot{\boldsymbol{r}_{2}}, \end{equation}

(5) \begin{equation}\boldsymbol{m}=\boldsymbol{Chaotic}\_\boldsymbol{value},\end{equation}

where ${\boldsymbol{r}_{1}}$ and ${\boldsymbol{r}_{2}}$ are random vectors in the range [0,1]. f is a convergence factor whose value decreases nonlinearly from 2.5 to 0 through the iteration process (in both the exploration and exploitation phases). $\boldsymbol{a}$ is a random vector that determines the distance between the chimps and prey in the range $\left[ -f,f \right]$ . $\boldsymbol{m}$ is a chaotic vector calculated based on various chaotic maps to represent the effect of the sexual motivation of chimps in the hunting process. $\boldsymbol{c}$ is the control coefficient used to describe the chimps’ expulsion and pursuit of prey, defined as a random vector in the range [0,2]. The process of chimps attacking the prey can be modelled using Eqs. (6)–(8).

(6) \begin{equation}\begin{aligned} &{{\boldsymbol{d}}_{Attacker}}=\left| {{\boldsymbol{c}}_{1}}{\boldsymbol{x}_{Attacker }}-{{\boldsymbol{m}}_{1}}\boldsymbol{x} \right|,{{\boldsymbol{d}}_{Barrier}}=\left| {{\boldsymbol{c}}_{2}}{\boldsymbol{x}_{Barrier}}-{{\boldsymbol{m}}_{2}}\boldsymbol{x} \right|\\&{{\boldsymbol{d}}_{Chaser}}=\left| {{\boldsymbol{c}}_{3}}{\boldsymbol{x}_{Chaser}}-{{\boldsymbol{m}}_{3}}\boldsymbol{x} \right|, {{\boldsymbol{d}}_{Driver}}=\left| {{\boldsymbol{c}}_{4}}{\boldsymbol{x}_{Driver}}-{{\boldsymbol{m}}_{4}}\boldsymbol{x} \right|,\end{aligned}\end{equation}

(7) \begin{equation}\begin{aligned} & {\boldsymbol{x}_{1}}={\boldsymbol{x}_{Attacker}}-{{\boldsymbol{a}}_{1}}({{\boldsymbol{d}}_{Attacker}}),\text{ }{\boldsymbol{x}_{2}}={\boldsymbol{x}_{Barrier}}-{{\boldsymbol{a}}_{2}}({{\boldsymbol{d}}_{Barrier}}) \\ & {\boldsymbol{x}_{3}}={\boldsymbol{x}_{Chaser}}-{{\boldsymbol{a}}_{3}}({{\boldsymbol{d}}_{Chaser}}),\text{ }{\boldsymbol{x}_{4}}={\boldsymbol{x}_{Driver}}-{{\boldsymbol{a}}_{4}}({{\boldsymbol{d}}_{Driver}}),\\\end{aligned}\end{equation}

(8) \begin{equation}\boldsymbol{x}(t+1)=\frac{{\boldsymbol{x}_{1}}+{\boldsymbol{x}_{2}}+{\boldsymbol{x}_{3}}+{\boldsymbol{x}_{4}}}{4}.\end{equation}

As indicated in Eqs. (6)–(8), after population initialisation, four optimal solutions are selected to define the positions of the attackers, barriers, chasers, and drivers. The positions of other chimps in the population are updated considering these positions. In other words, the position of the prey is estimated by those of the four best individuals, and the locations of other chimps are randomly updated in the vicinity.

In the last stage of hunting, the individuals achieve food satisfaction and exhibit natural wild behaviour in the presence of social opportunities. Specifically, the chimps attempt to obtain food by force and chaos. The chaotic behaviour of chimps in the final stage helps alleviate the problems of local optimum trapping and low convergence speed, and chaotic maps can be used to enhance the ChOA performance. The updating process is modelled using Eq. (9), where $\mu$ is a random number in [0,1]. Because certain chimps may not have any sexual motivation in the hunting process, a probability of 50% is set to determine whether the position update strategy of the chimps is normal (Eq. (2)) or chaotic (Eq. (9)). The process flow of the ChOA is illustrated in Fig. 1.

(9) \begin{equation}{\boldsymbol{x}_{chimp}}(t+1)=\left\{ \begin{aligned}& {\boldsymbol{x}_{\text{prey}}}(t)-\boldsymbol{a}\cdot\boldsymbol{d} \qquad if \ \mu< \text{0}\text{.5} \\[5pt] & \boldsymbol{Chaotic}\_{\boldsymbol{value}}\qquad if \ \mu> \text{0}\text{.5} \\\end{aligned} \right..\end{equation}

Figure 1. Process flow of ChOA.

2.2 Improved ChOA

Although the ChOA has many advantages and has been applied in many fields, it has a few limitations. Specifically, the population initialisation is based on random distribution, which leads to a low population diversity, aimlessness in the individual early search, and low convergence speed of the algorithm. In addition, because the positions of the other chimps are updated based on those of the attackers, these chimps may be trapped in the local optimum, which limits their exploration of new areas in the search space owing to their solution space being concentrated around the attackers’ solutions.

Consequently, the original approach is not suitable for solving highly complex optimisation problems owing to limitations such as low diversity, vulnerability to local optimum trapping, and low convergence speed. To address these limitations, the population initialisation and location update methods are enhanced in this study. To enhance the exploitation phase of the standard ChOA, the search process equations are modified, as described in the following subsections.

2.2.1 Sobol sequence for population initialisation

For all algorithms involving sampling, a superior sample distribution can be obtained using evenly distributed random numbers. In the original ChOA, the population initialisation process is based on random distribution, which leads to a low diversity and aimlessness of individual optimisation. In the improved algorithm, the Sobol sequence is used to initialise the population. The dimensions of the Sobol sequence are based on radical inversion with a factor of 2; however, each dimension has a specific generation matrix. Notably, the Sobol sequence is a low discrepancy sequence, in which a deterministic pseudo-random number sequence is used to replace the pseudo-random number sequence to fill the sample space as evenly as possible. Therefore, the Sobol sequence is associated with a high computational efficiency and wide coverage of sampling points in solving the initial sampling problems.

Considering the optimisation problem of interest, the upper and lower limits of the two-dimensional space are set as $\boldsymbol{lb}=[\text{5 25,000}]$ and $\boldsymbol{ub}=[\text{6 30,000}]$ to describe the cruise state. The random number in the sampling space is generated using Eq. (10), and the population size is set as 30. ${S}_{i}$ is the ith random number generated by the Sobol sequence in the range [0,1].

(10) \begin{equation}{{\boldsymbol{X}}_\boldsymbol{n}}=\boldsymbol{lb}+{{S}_{i}}\cdot (\boldsymbol{ub-lb}).\end{equation}

2.2.2 Weight factor

To alleviate the problem of the other chimps falling into the local optimum, the iwChOA incorporates a position-weighted relationship based on a weight factor that can be mathematically formulated as follows:

(11) \begin{equation}{{\beta }_{1}}=\frac{{{\left| \boldsymbol{x} \right|}_{1}}}{{{\left| \boldsymbol{x} \right|}_{1}}+{{\left| \boldsymbol{x} \right|}_{2}}+{{\left| \boldsymbol{x} \right|}_{3}}+{{\left| \boldsymbol{x} \right|}_{4}}}, \end{equation}

(12) \begin{equation}{{\beta }_{2}}=\frac{{{\left| \boldsymbol{x} \right|}_{2}}}{{{\left| \boldsymbol{x} \right|}_{1}}+{{\left| \boldsymbol{x} \right|}_{2}}+{{\left| \boldsymbol{x} \right|}_{3}}+{{\left| \boldsymbol{x} \right|}_{4}}}, \end{equation}

(13) \begin{equation}{{\beta }_{3}}=\frac{{{\left| \boldsymbol{x} \right|}_{3}}}{{{\left| \boldsymbol{x} \right|}_{1}}+{{\left| \boldsymbol{x} \right|}_{2}}+{{\left| \boldsymbol{x} \right|}_{3}}+{{\left| \boldsymbol{x} \right|}_{4}}},\end{equation}

(14) \begin{equation}{{\beta }_{4}}=\frac{{{\left| \boldsymbol{x} \right|}_{4}}}{{{\left| \boldsymbol{x} \right|}_{1}}+{{\left| \boldsymbol{x} \right|}_{2}}+{{\left| \boldsymbol{x} \right|}_{3}}+{{\left| \boldsymbol{x} \right|}_{4}}},\end{equation}

(15) \begin{equation}\boldsymbol{x}(t+1)=\frac{1}{{{\beta }_{1}}+{{\beta }_{2}}+{{\beta }_{3}}+{{\beta }_{4}}}\times \frac{{{\beta }_{1}}{{\left| \boldsymbol{x} \right|}_{1}}+{{\beta }_{2}}{{\left| \boldsymbol{x} \right|}_{2}}+{{\beta }_{3}}{{\left| \boldsymbol{x} \right|}_{3}}+{{\beta }_{4}}{{\left| \boldsymbol{x} \right|}_{4}}}{4}.\end{equation}

The learning rate is introduced in Eq. (8) to obtain Eq. (15). ${{\beta }_{1}}$ , ${{\beta }_{2}}$ , ${{\beta }_{3}}$ , and ${{\beta }_{4}}$ are the learning rates of other chimps with respect to the attacker, barrier, chaser and driver, respectively, in terms of the Euclidean distances. The dynamic learning rate in Eq. (15) varies across iterations, which helps increase the convergence speed and ensure the avoidance of the local optima in scenarios in which the attackers, barriers, chasers and drivers are less likely to know the position of the prey.

3.1 Dynamics model

For simplicity, the following assumptions are implemented:

1. (1) The curvature of the earth is ignored, and the ground coordinate system is the inertial coordinate system.

2. (2) The instantaneous equilibrium hypothesis is applied.

3. (3) Only longitudinal motion is considered.

Equations (16)–(21) represent a typical dynamic model of a hypersonic vehicle:

(16) \begin{equation}\frac{dV}{dt}=\frac{T\cos \alpha -D}{m}-g\sin \gamma, \ {\rm }\end{equation}

(17) \begin{equation}\frac{d\gamma }{dt}=\frac{L+T\sin \alpha }{m}-\frac{g}{V}\cos \gamma, \ {\rm }\end{equation}

(18) \begin{equation}\frac{d\alpha }{dt}=q_d-\frac{d\gamma }{dt}, \ {\rm }\end{equation}

(19) \begin{equation}\frac{dx}{dt}=V\cos \gamma, \ {\rm }\end{equation}

(20) \begin{equation}\frac{dy}{dt}=V\sin \gamma, \ {\rm }\end{equation}

(21) \begin{equation}\frac{dm}{dt}=-T/(g\cdot isp),\end{equation}

where $g={{g}_{0}}{{(\dfrac{{{R}_{e}}}{r})}^{2}}$ is the acceleration due to gravity, ${{g}_{0}}$ is the acceleration due to gravity on the ground, and ${{R}_{e}}$ is the radius of the Earth. D and L are the drag and lift forces, respectively, calculated as

(22) \begin{equation}\\D={{q}_{d}}{{S}_{r}}{{C}_{D}},\end{equation}

(23) \begin{equation}\\L={{q}_{d}}{{S}_{r}}{{C}_{L}}.\end{equation}

The dynamic pressure ${{q}_{d}}$ can be defined as

(24) \begin{equation}\\{{q}_{d}}=\frac{1}{2}\rho {{v}^{2}},\end{equation}

$\rho$ is the atmospheric density, which can be calculated as

(25) \begin{equation}\\\rho ={{\rho }_{0}}{{e}^{((r-{{R}_{e}})/{{h}_{0}})}},\end{equation}

where ${{\rho }_{0}}$ is the atmospheric density at sea level (1.225570827014494), and ${{h}_{0}}$ is the elevation at which the density is being evaluated (7254.24).

3.2 Hypersonic cruise vehicle model

The considered model adopts three engines in hypersonic cruise flight, which operate in a combined cycle propulsion system. Such power plants can cover the mission requirements in different speed ranges such as low speed, supersonic and hypersonic regimes [Reference Keshmiri and Colgren33]. The ratio of the engine thrust to the Mach number can be calculated using the following analytical formulas.

(26) \begin{equation}\left\{\begin{aligned} 0.0\le &\;Mach\le 2.0 \\ T=&\;PLA*(2.99*{{10}^{5}}-32.8*(h)+1.43*{{10}^{-3}}*{{(h)}^{2}}-2.29*{{10}^{-8}}*{{(h)}^{3}}\\ &\; +3.75*{{10}^{3}}*{{(Mach)}^{3}}) \\\end{aligned},\right.\end{equation}

(27) \begin{equation}\left\{\begin{aligned} 2.0\le& \;Mach\le 5.0 \\ T=&\;PLA*(396.103*Mach-703.1*Mac{{h}^{2}}+816.19*Mac{{h}^{3}}-442.48*Mac{{h}^{4}} \\ &+118.28*Mac{{h}^{5}}-15.34*Mac{{h}^{6}}+0.772*Mac{{h}^{7}})*5*{{10}^{3}};\end{aligned},\right.\end{equation}

(28) \begin{equation}\left\{\begin{aligned} 5.0\le&\; Mach \le 6.0 \\T=&\;PLA*(396.103*Mach-703.1*Mac{{h}^{2}}+816.19*Mac{{h}^{3}}-442.48*Mac{{h}^{4}} \\ &+118.28*Mac{{h}^{5}}-15.34*Mac{{h}^{6}}+0.772*Mac{{h}^{7}})*5*{{10}^{3}};\end{aligned}.\right.\end{equation}

The change in the hypersonic vehicle mass during flight, which depends on the thrust, gravitational acceleration and specific impulse, can be obtained using Eq. (21). The specific impulse can be obtained via interpolation of the data in Table 1.

Table 1. Engine specific impulse values

The aerodynamics parameters are obtained via piecewise calculation [Reference Keshmiri and Colgren33]:

(29) \begin{eqnarray}{{C}_{D}}={{f}_{D}}(Mach,\alpha ), \end{eqnarray}

(30) \begin{eqnarray}{{C}_{L}}={{f}_{L}}(Mach,\alpha ), \end{eqnarray}

where ${{f}_{D}}$ and ${{f}_{L}}$ are piecewise interpolation functions to analyse the drag and lift coefficients, respectively.

4.1 Objective function

As a key phase of the trajectory, the cruise phase has its own objectives. The take-off performance optimisation problem of a hypersonic vehicle pertains to determining the optimal cruise speed and altitude, ${\boldsymbol{S}^{*}}(t)=[{{M}^{*}}(t),{{A}^{*}}(t)]$ , to minimise the take-off mass subject to range requirements and other constraints. Typically, the objective function is required to be maximised. Therefore, the negative take-off mass is set as the performance index:

(31) \begin{equation}J(s)=-\left| {{mass}_{1}}+{{mass}_{2}} \right|,\end{equation}

where ${{mass}_{1}}$ is the weight of the fuel consumed in the climb phase, and ${{mass}_{2}}$ is the initial mass of the vehicle in the cruise phase. The optimisation objective is to determine the best cruise state to maximise the objective function specified in Eq. (31), subject to the 3-DOF dynamics presented in Eqs. (16)–(21) and constraint conditions.

4.2 Constraints

The trajectory optimisation of a hypersonic vehicle in the climb phase is a typical nonlinear optimal control problem with state and control constraints. The Gaussian pseudo-spectral method can obtain a superior solution with a lower computational cost. Therefore, the flight process in the climb phase is modelled using the Gauss pseudo-spectral method [Reference Williams, Sun, Bai and Yan34]. In the optimisation process, the initial speed and altitude in the climb phase are set as ${{v}_{0}}=100$ m/s and ${{h}_{0}}=50$ m, respectively. The end of the climb phase marks the start of the cruise phase, which is the search space in this study. In this phase, the Mach number is $M=[5,6]$ and altitude is $A=[25,000,30,000]$ . The search space is divided equally, and certain points are selected as the end points of the climb phase. The Gauss pseudo-spectral method is used to determine the changes in the flight state variables in the climb process, with the variable of interest being the vehicle mass.

As mentioned earlier, the objective of the considered problem is to determine the optimal cruise speed and altitude to minimise the aircraft take-off mass. The main requirement for aircraft control in the cruise phase is to maintain stable flight and satisfy the range requirements. Because the state variables in the cruise phase change gradually, flight control can be realised in a facile manner. The stability of the trajectory in the cruise phase must be ensured through guidance control, which mainly involves altitude and speed control, the control variables of which are the angle-of-attack and throttle equivalent, respectively. In this study, the classic proportional integral–derivative (PID) control is adopted for cruise control, although the details are not presented in this paper.

The cruise flight process contains four control constraints and two path constraints, as indicated in Eq. (32). All constraints except for one are inequality constraints. The equality constraint ensures that the vehicle attains the cruise range. The inequality constraints ensure that the change in the angle-of-attack and throttle equivalent do not exceed the pre-programmed boundaries, and the rate of change of the attack angle does not exceed the limits. The range, angle-of-attack and rate of change of the angle-of-attack are expressed in units of m, ${\deg}$ and deg/s, respectively, and the fuel equivalent is expressed as a percent point. The path constraint ensures that the dynamic pressure and heat flux density of the vehicle do not exceed the maximum limits.

(32) \begin{equation}\left\{ \begin{matrix} { }Range=1.3*{{10}^{7}}\hfill \\ { }0\le \alpha \le 10\hfill,-3\le \dot{\alpha }\le 3,0\le PLA\le 100\%\hfill \\ { }q=\frac{1}{2}\rho {{V}^{2}}\le {{q}_{\max }},{{{\dot{q}}}_{s}}={{k}_{p}}{{\rho }^{0.5}}{{V}^{3.08}}\le {{{\dot{q}}}_{s\max }} \\\end{matrix} \right..\end{equation}

4.3 Optimisation based on iwChOA

The take-off performance optimisation process based on the iwChOA is illustrated in Fig. 3. First, the parameters of the aircraft and proposed algorithm are initialised, and the fitness of the initial population is calculated. The iwChOA continuously updates the position according to the process described in Section 2.2.3 until the optimisation result satisfies the iteration termination conditions. The fitness value is obtained by calculating the flight state of the aircraft, which is the take-off mass when the flight range is attained.

Figure 3. Optimisation of the take-off mass.