An Investigation of the Bifurcation Behavior of an F-18 Aircraft Model

Bifurcations of equilibria and limit cycles of an F-18 aircraft model have been investigated in this paper based on one or two continuation parameters. First, it is shown that the aircraft can be in a stable straight-and-level flight state by coordinating the elevator deflection and the engine thrust and fixing the other parameters. Second, the aircraft may be disturbed by perturbations out of the longitudinal plane and get into the lateral-directional motion mode by a branch point bifurcation. Third, when the straight flight state loses its stability, body axis angular rates will present periodic or quasi-periodic motion pattern by the Hopf or Neimark–Sacker bifurcation. Finally, bifurcation structure of limit cycles has been exhibited, including the generalized Hopf bifurcation, the Neimark–Sacker bifurcation, the Hopf-Hopf bifurcation and the fold-Neimark–Sacker bifurcation. Meanwhile corresponding bifurcation curves in two-parameter plane have been depicted with the help of numerical continuation techniques.


Introduction
As a modern means, aircrafts are playing an increasingly important role in transportation and combat industry. During the flight, it is necessary to maintain the stability, controllability and safety of the aircrafts to complete transportation and combat tasks quickly and efficiently. Before the start of the mission, the prediction of aircraft trajectories provides important information for the risk assessment of air transportation safety. Accurate trajectory prediction can be achieved through the actual combat drill in the early stage. However, due to existence of the large number of aircrafts in the airspace, the drilling simulation costs can be expensive or even restricted. Thus, a computationally efficient method is needed for the simulation of aircraft dynamics in order to achieve the real time prognostics of the air transportation system [1]. In order to elucidate the dynamics of the aircrafts, the usual approach is to develop the kinematics and dynamics equations of the barycenter of aircrafts based on a certain coordinate system. A high fidelity simulation model provides engineers with the numerical power to test new airplane designs or any modifications to existing ones in a simulation environment prior to flight tests [2].
The kinematic equations of the aircrafts involve high nonlinearity because of the strong coupling between the rigid-body and structural dynamics of the aircrafts, including structural, aerodynamic, or control nonlinearities. Nonlinear problems in aircraft flight dynamics and control have been well recognized and widely documented ever since the dawn of aviation [3]. Bifurcation and continuation methods have been identified as a practical tool to analyze the dynamics of the aircrafts in presence of nonlinearities, which can be used to discern the bifurcation points, construct the bifurcation diagrams, describe the parameter plane structure diagrams, and so on [4][5][6][7][8][9].
In [10], the authors study the pitch attitude dynamics of an asymmetric magnetic spacecraft in an almost circular orbit under the influence of a gravity gradient torque. It is shown that the pitch motion with perturbations exhibits heteroclinic chaotic behavior by means of the Melnikov method. Some periodic pitch motions persist in the perturbed system with the same period as the orbital motion of the spacecraft. In [11], the authors evaluate the effects of flexibility on the dynamics, stability and control of elastic aircraft using an analysis framework based on bifurcation and continuation tools, which reveals the nonlinear dynamic coupling between rigid-body and flexible modes in regions of the flight envelope. Furthermore, the efficiency of the approach relative to traditional nonlinear simulation is discussed and the robustness of idealized control law designs to unmodeled elastic modes is also investigated. In [12], the author explores the impact of kinematic nonlinearities on the dynamics of a highly deformable cantilevered wing. The impact of the degree of freedom on the stability of nonzero trim states has been examined and limit-cycle oscillations are also investigated in terms of amplitude of vibration. The effect of variation of system parameters such as stiffness ratio, aspect ratio and root angle of attack has also been studied.
Based on bifurcation analysis method, the problem of deep-stall recovery has been studied for the aircraft without longitudinal static stability in [13], where a finite-time prescribed performance deep-stall recovery control scheme is designed and stability of the closed-loop system is proved by the common Lyapunov functional method. Rohith and Sinha in [14] examine possible routes to chaos in the post-stall dynamics of an F-18 high-model with the wind as the driving agent. It has been established that there are three routes to chaos including quasi-periodic, period-doubling and intermittency route according to different post-stall flight conditions. The flight performance and dynamics analysis of F-18 aircraft model have been investigated in [15,16]. In [15], the problem of how to assess the safe landing capability after encountering impairment has been studied, and performance parameters and stability along the manoeuvre performance curves are obtained by the constrained bifurcation and continuation analysis techniques, which detailedly reveals how flight performance and stability change due to the aircraft impairment. In [16], the Hopf bifurcation and the pitchfork bifurcation in longitudinal flight state are discussed as well as the effect of aileron deflections on the turn properties are also studied. For longitudinal flight, the existence of periodic orbits and equilibria at a nonzero bank angle present a risk for the flight safety in a situation. For the aircraft turns, the fold bifurcations and associated jumps reveal the effective range of the lateral control.
Even though many studies about the F-18 aircraft model have come forth, however, bifurcations of equilibria and limit cycles about this model have not been adequately researched, such as the cusp and the zero-Hopf bifurcation of equilibria, bifurcation structure of limit cycles in two-parameter plane and detection of strong and weak resonance points. As a further research, several novel issues have been investigated in this paper. First, the transition mechanism of equilibrium states from longitudinal flight to lateral-directional flight have been elaborated. Second, relatively comprehensive bifurcations of limit cycles have been investigated and corresponding bifurcation curves have been depicted in two-parameter planes for two pairs of parameters, where symmetrically distributed structure have been presented. Third, strong and weak resonance points and fold-Neimark-Sacker points have been detected on the Neimark-Sacker bifurcation curve of limit cycles, near which multiple period limit cycles and two-dimensional tori can emerge. In this paper, bifurcation analysis and continuation methods are implemented in the MATLAB environment using the dynamical system software package Matcont [17].
This paper is organized as follows. An aerodynamic model of the F-18 aircraft is introduced in Sect. 2. The equilibria in longitudinal motion and lateral-directional motion as well as possibly associated bifurcations have been expounded in Sects. 3 and 4, respectively. The bifurcation structure of equilibria and limit cycles in twoparameter planes has been presented in Sect. 5 and corresponding phase diagrams have also been displayed. The conclusions are given in Sect. 6.

The Kinematic Equation for the Aerodynamic Model F-18
To describe aircraft flight dynamics, three common coordinate systems are introduced here. One is the inertial reference frame O E -X E Y E Z E , and the origin O E is located at some given fixed point on the ground. The axis O E X E is in a horizontal plane and points to some fixed direction; the axis O E Z E is vertically downward; the axis O E Y E is perpendicular to the plane O E X E Z E , following the right-hand rule. Another is the body- Here the aircraft is assumed to be an absolutely rigid body with the constant mass and a symmetric plane during the flight, and the earth is flat with no rotating. The time rate of change of mass and inertia is assumed negligible, and fuel-sloshing effects are ignored. Assuming that the body-fixed axis system coincides with the principal axes of the aircraft, and the engine alignment and thrust vector through the centre of gravity are along the X axis of the aircraft. The effect of the engine thrust and three types of control deflections-the elevator, the aileron and the rudder have been considered, while the effect of other types of flight control surfaces, such as flaps, slats, spoilers and stabilizers, has been ignored. Some more assumptions regarding the geometry of the aircraft can refer to the reference [9].
Based on the above coordinate systems and assumptions, the dynamical model of the F-18 aircraft has been established as the following nonlinear differential equations [9,21].
(2.1g) = p + qtan sin + rtan cos , The Mach number M a is defined as the ratio of the aerodynamic velocity V to the sound velocity a at sea level. There are two aerodynamic angles between the body and the velocity coordinate system, the angle of attack and the sideslip angle , which are shown in Fig. 1a. Here is restricted to a low-angle-of attack range about between −5 and 40 degrees, which can guarantee that the lift is greater than the drag; the relatively low sweep wing works well; the model is symmetric with respect to the longitudinal plane. Three Euler angles between the body and the inertial reference frame, the bank angle , the pitch angle and the body yaw angle , which are shown in Fig. 1b. Three body axis angular rates-roll p, pitch q and yaw r are displayed in Fig. 1c. Two flight-path angles between the velocity and the inertial reference frame are defined here, the wind-axis roll angle and the angle of inclination , which are exhibited in Fig. 1d. In fact, there is also a flight-path azimuth angle , which has not been used or defined here. Equations (2.1a-2.1c) are the translational dynamics and (2.1d-2.1f) are the rotational dynamics and the last three are the rotational kinematics. T is the engine thrust and the throttle = T∕T max is the fraction of the maximum available engine thrust. The constant moments of inertia I x , I y and I z are listed in Tab.1 in Appendix. The sin = cos cos sin − sin sin cos − sin cos cos cos , sin cos = sin cos sin + sin cos cos − sin sin cos cos , cos cos = sin sin + cos cos cos . Fig. 1 a Aerodynamic angles-and . b Euler angles-, and . c Angular rates-p, q and r. d Flight-path anglesand geometric datum and aerodynamic parameters C D , C Y , C L , C m , C l and C n are listed in Tab.2 in Appendix, which are presented in first or second order Taylor series forms based on the symmetrical aircraft configurations [9,20,21]. About this model, the nonlinear factors stem from the kinematic coupling, inertia coupling, nonlinear aerodynamics and so on, which result in nonlinear dynamics, such as loss of stability, stall and post-stall [9].

The Equilibria in Longitudinal Motion
As the yaw angle appears in constrained Eq. (2.1i) and has no effect on the other variables, the Eqs. (2.1a-2.1h) are only considered here. An aircraft usually flies in the conditions that are constrained such that more than one control is required to be employed. Usually, the flight states can be uncoupled into longitudinal motion and lateral-directional motion. At first, the longitudinal trim and its stability will be analyzed here.

Straight and Level Flight Trim
In this flight state, the constraints of aircraft equations along this motion are given as     right column of Fig. 2. The left column of Fig. 2 indicates that, with the increase of , the existence of straight and level flight states demands that the throttle and the upward elevator deflection e are gradually increasing while the Mach number M a is gradually decreasing.

Climbing Flight in the Vertical Plane
From Fig. 2a-c, it is shown that, for a given , there exists a unique point ( e , M a ) such that the aircraft is in stable level flight in longitudinal plane. When e varies with the other parameters fixed, the level flight state will be violated. Figure 3a And the magnitude of increase in is larger than the one in , which indicates that the aircraft from nose down to nose up is climbing in vertical plane with a deceleration. By analysis of eigenvalues of the equilibrium curve in Fig. 3a-c, it is shown that when e is chosen in the right of H 2 ( e > −0.1932rad), the equilibria are stable; when e is chosen between H 1 and H 2 (−0.2438rad < e < −0.1932rad), the equilibria are unstable; when e is chosen in the left of the point H 1 ( e < −0.2438 rad), the equilibria are stable. Two typical phase diagrams for e =−0.17rad and e =−0.2rad are exhibited in Fig. 4a, b, respectively.  Fig. 3d. A red time-series diagram for = 0.9, a = 0.003662rad, r = −0.002937rad and e = −0.2rad, which is corresponding to the red branch in Fig. 3d 1 3

The Equilibria in the Lateral-Directional Dynamic Modes
In contrast with the fold bifurcation, the transcritical bifurcation and the pitchfork bifurcation are regarded as the branch point bifurcation [22]. Figure 3a-c show that with the variation of e , two branch points BP 1 and BP 2 will be encountered, at which the Jacobian matrix M of the system (2.1a-2.1h) has only one zero eigenvalue and from which other equilibrium branches will emerge. In order to discriminate the direction of the secondary branch to switch branches, the left and the right unit eigenvectors of M need to be computed [22], and the original free parameter e becomes the branch parameter as well as two other parameters need to be activated [18]. Here e , a and r are chosen as the active parameters to obtain the other equilibrium branches. The bifurcation structures near the points BP 1 and BP 2 are exhibited in Fig. 3d. From the point BP 2 , besides the longitudinally stable equilibrium state (blue branch), two other stable lateral-directional equilibrium branches (red and magenta) emerge, which indicates that BP 2 is a pitchfork bifurcation point, from which three stable equilibrium branches bifurcate. The red branch will be obtained by choosing e as a branch parameter and e , a and r as active parameters. However, the magenta branch will be acquired, if e is still chosen as the continuation parameter with = 0.9, a = 0 , r = 0 and the initial values are different from those in the blue branch. This shows that there is a longitudinal equilibrium state and a lateraldirectional equilibrium state coexisting for = 0.9, a = 0 and r = 0 as well as e ∈ (−0.1932, −0.1562) as a continuation parameter.
If e , a and r are all activated, three equilibrium states will be coexistent. Three typical time-series diagrams for e = −0.18rad from different initial values have been plotted in Fig. 4c. For e in the left of H 2 , the blue branch loses its stability and the stability of the other two branches remain unchanged. Two typical time-series diagrams for e = −0.2rad have been displayed in Fig. 4d.
H 1 and H 2 in Fig. 3a-c are subcritical Hopf bifurcation points, whose location has been detected by the test function H (M a , , ; e ) = ∏ i>j ( i (M a , , ; e ) + j (M a , , ; e ))) in Matcont. Wherein i and j are the eigenvalues of the Jacobian matrix M of the system (2.1a-2.1h). The function H is real and smooth and has a regular zero at H 1 or H 2 with 1,2 = ±i 0 . In the numerical analysis of bifurcation, the bialternate product (2M ⊙ I n ) has been introduced to compute H and the test function can be expressed as H =det 2M ⊙ I n . In addition, according to the Hopf bifurcation theory in [23], the first Lyapunov coefficient can be obtained by deriving the normal form of Hopf bifurcation. It is shown that the first Lyapunov coefficients are positive for H 1 and H 2 . From them, two series of unstable limit cycles bifurcate towards the opposite directions with the variation of e , which is shown in Fig. 5a. The subcritical Hopf bifurcations may show jump phenomena to either a steady state or another (stable) limit cycle branch. At the point BP 1 , three parameters require to be activated to get the additional equilibrium branches. When e , a and r are chosen as active parameters and e as the branch parameter, a stable black branch bifurcates from the point BP 1 and another branch point BP 3 will be encountered. Next, if e varies and = 0.9, a = 0.00313rad and r = −0.0032rad, an unstable green branch and a purple branch will be emanated from the point BP 3 with different initial values. On the green branch, a subcritical Hopf bifurcation point H 3 will be detected, from which a series of unstable limit cycles shown in Fig. 5c bifurcate towards the right with the variation of e . At e = −0.1927 rad, a Neimark-Sacker bifurcation point NS will be encountered, where the Poincaré map associated with the periodic orbit has a simple pair of conjugate complex Floquet multipliers with moduli 1. Generically this corresponds to a bifurcation of an invariant torus, on which the flow of the system presents periodic or quasi-periodic motion. When the NS point is used as a starting point for the computation of a discrete orbit in the Poincaré section for a slight variation of e , then after a transient the timeseries exhibits modulated oscillations with two frequencies near the original limit cycle. A typically stable two-torus is displayed in Fig. 5d. In addition, the purple branch is divided into a stable branch (upper) and an unstable one (lower) by the limit point LP 1 . For = 0.9 , a = 0.00313rad, r = −0.0032 rad and e = −0.20769rad, two coexistent equilibrium states (the upper purple branch and the black branch in Fig. 3d) from different initial values are demonstrated in Fig. 5b.

Bifurcation Curves in Two-Parameter Plane
In previous sections, stability and bifurcation of the equilibria in one-parameter plane have been discussed. In this section, bifurcation structure of equilibria and limit cycles in two-parameter plane will be analyzed. Figure 5a is a single parameter Hopf bifurcation diagram for = 0.9 , a = 0 and r = 0 . In order to obtain Hopf bifurcation curves and deuterogenic bifurcations, two-parameter bifurcation needs to be discussed.

Bifurcation Structure in ı e -ı a Plane
For r = 0 , the Hopf bifurcation curve from H 1 is depicted in Fig. 6 with a and e as continuation parameters. It is shown that as the cyan curve l 1 extends upwards and downwards from the point H 1 , symmetrically generalized Hopf points GH i (i = 1, ⋅ ⋅ ⋅, 8) and Hopf-Hopf HH j ( j = 1, ⋅ ⋅ ⋅, 4) can be detected. At each point GH i , the first Lyapunov coefficient is zero and a fold bifurcation curve of limit cycles will bifurcate from it. At each point HH j , there are two pairs of purely imaginary eigenvalues in the Jacobian matrix M and a Neimark-Sacker bifurcation curve Fig. 6 Bifurcation structure diagram in e − a plane for = 0.9 and r = 0 of limit cycles will emanate from it. Owing to the similarity, only the bifurcation structure near the points HH 2 and GH 3 have been demonstrated as representatives.
It is shown that there is a Neimark-Sacker bifurcation curve l 2 bifurcating from the point HH 2 , which happens to end at the point HH 3 . On the curve l 2 , two pairs of symmetrically strong resonant points R 2 and R ′ 2 as well as R 4 and R ′ 4 will be monitored, where a period-2 and a period-4 limit cycle will emerge, respectively [19]. Two representative phase diagrams for the points R 2 and R 4 have severally been exhibited in Fig. 7a, b. In addition, a pair of weak resonant points R 9 and R ′ 9 as well as two symmetrical Fold-Neimark-Sacker bifurcation points FNS and FNS ′ can also be detected. Therein a period-9 limit cycle and a two-dimensional torus will appear, which have been displayed in Fig. 7c, d, respectively. The appearance of multi-periodic closed orbits and the torus shows that the pitch angle, the yaw angle and the roll angle vary in periodic, multi-periodic or quasi-periodic pattern.
For the point GH 3 , from which a fold bifurcation curve of limit cycles l 3 bifurcates. When the parameter pair ( e , a ) has been chosen above the curve l 3 , a stable (a) (b) (c) (d) Fig. 7 The phase diagrams for the points in Fig. 6. a A period-2 limit cycle for the point R 2 with ( e , a ) = (−0.261254rad, 0.022786rad) . b A period-4 limit cycle for the point R 4 with ( e , a ) = (−0.256185rad, 0.010204rad) . c A two-dimensional torus for the point FNS with ( e , a ) = (−0.263763rad, 0.025402rad) . d A period-9 limit cycle for the point R 9 with ( e , a ) = (−0.254568rad, 0.131558rad) equilibrium and a stable limit cycle coexist with an unstable limit cycle between them. When the parameter pair ( e , a ) crosses the curve l 3 , a fold bifurcation of limit cycles occurs and two limit cycles collide to disappear with the stable equilibrium left. The phase diagrams for the points A and B have been displayed in Fig. 8a, b, respectively.

Bifurcation Structure in ı e -ı r Plane
For a = 0 , the Hopf bifurcation curve from the point H 1 has been shown in Fig. 9 with e and r as continuation parameters, which forms a symmetric figure-eight pattern and terminates at the other Hopf bifurcation point H 2 by coincidence. With  the extension of the Hopf curve, two zero-Hopf points ZH 1,2 will be encountered, at which the system has an equilibrium with one zero eigenvalue and a pair of purely imaginary eigenvalues ±i 0 . The zero-Hopf bifurcation here corresponds to the case for s = −1 and < 0 [23]. From the points ZH 1,2 , two fold bifurcation curves of equilibria will emerge and converge at the point CP to form a small wedge region. It is remarkable that here the point CP just happens to be the point BP 1 in Fig. 3. Inside the small wedge, there is only a stable equilibrium, which is shown in Fig. 10b for the parameter pair P � (−0.2037rad, 0) . The point P ′ corresponds to = 0.9, a = 0 , r = 0 and e = −0.2037rad , with which there are two branches (magenta and blue) associated in Fig. 3d. When e = −0.2037rad , only the upper magenta branch is stable, which happens to be consistent with the result in Fig. 10b. Outside the small wedge, there are three equilibria, two stable and one unstable. Two stable equilibria C and D for the parameter pair P (−0.2037rad, 0.0001rad) have been exhibited in Fig. 10a, where the unstable equilibrium has not been presented.

Conclusions
Bifurcations of equilibria and limit cycles of an F-18 aircraft model have been investigated in this paper. It is shown that the aircraft can deviate from the longitudinal plane by the branch point bifurcation and then enter into the lateral-directional straight motor pattern with the swing of the rudder and ailerons. When the equilibria lose their stability, the aircraft can enter into periodic or quasi-periodic motion mode by Hopf bifurcation or Neimark-Sacker bifurcation, respectively. In ea parameter plane, the strong and weak resonant points can be detected on the Neimark-Sacker (a) (b) Fig. 10 a Two stable equilibria C and D for the parameter pair P (−0.2037rad, 0.0001rad) in Fig. 9 with different initial values. b A stable equilibrium E for the parameter pair P � (−0.2037rad, 0) in Fig. 9 bifurcation curve, near which period-2, period-4 and period-9 trajectories can emerge. In er parameter plane, it is revealed that two stable equilibria coexist and then overlap to form one equilibrium through the cusp bifurcation. These results can provide a theoretical basis for the F-18 aircraft to expound and predict the transition mechanism of dynamics. Based on these, the elevator, the aileron and the rudder deflection angle can be coordinated such that the aircraft enters into the expected states.
There are some remarks about the treatments and results of this paper. (1) Generally, a high-dimensional system with nonlinearity can be reduced to a lowerdimensional one by some order reduction methods, as has been demonstrated in literature [24][25][26][27], where the system which needs to be reduced has a unified form MZ + CŻ + KZ = F . According to the feature of the model (2.1a-2.1i) in this paper, it is difficult or almost impossible to obtain the linear system and corresponding eigenvectors due to the aerodynamic parameters with multi-ply piecewise functions. So the generically theoretical analysis in the order reduction methods has not been considered and discussed here. (2) The treatments of this paper are carried out like this. First, some constraints have been imposed and the implicit expressions of the equilibria in longitudinal plane have been deduced. Next, the other equilibrium branches in lateral-directional plane are obtained by numerical simulation continuation methods. Finally, the relatively complete bifurcation structures in two-parameter plane have been presented in Figs. 6 and 9 with the help of the bifurcation analysis toolbox.  Table 2 The expressions of aerodynamic parameters for F-18 model Availability of data and material Not applicable.
Code availability Not applicable.

Conflicts of interest
The authors declare that there are no conflicts of interests regarding the publication of this manuscript.
Ethics approval Not applicable.
Consent to participate Not applicable.

Consent for publication Not applicable.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.