Abstract

The development of control systems for aerial vehicles necessitates a meticulous examination of their dynamic behavior. This research delves into an in-depth investigation of the dynamic behavior of the F-16 aircraft, employing refined mathematical models to analyze both its longitudinal and lateral motions, as well as their corresponding modes. These mathematical models are formulated in two conventional representations: state space equations and transfer functions. By utilizing these mathematical representations, two displacement autopilots have been developed, consisting of a pitch attitude autopilot based on the longitudinal equations and a roll attitude autopilot designed using the lateral equations. Proportional Integral Derivative (PID) controllers, encompassing inner loops, as well as Linear Quadratic Controllers (LQR), have been recruited as control system units. The control structures have undergone analysis utilizing Simulink models. The analyses have yielded favorable damping characteristics and faster responses in both longitudinal and lateral movements and modes.

Keywords

• Flight Dynamics
• Longitudinal Motion
• Lateral Motion
• Displacement Autopilots
• LQR
• PID
• F-16 Aircraft

1. Introduction

The F-16 aircraft represents a multi-role fighter renowned for its exceptional maneuverability. Its prowess has been demonstrated in both air-to-air and air-to-surface missions. This high level of agility is attained by intentionally shifting the center of gravity(CG)  from the stable region to the unstable region, resulting in a state of relaxed static stability (Reichert, 1993). Whenever the CG position falls within the unstable region, the F-16 aircraft can only sustain flight by relying on a flight control system known as the Control Augmentation System (CAS). This sophisticated flight control system alters the dynamic characteristics of the F-16, ensuring its stability even when the CG position is situated within the unstable CG region. Moreover, this control system offers the pilot the advantage of choosing task-specific control laws. For instance, a specialized control augmentation system is indispensable for high-performance fighter aircraft like the F-16, as it allows the pilot to execute intricate maneuvers, pushing the aircraft to its performance limits while performing tasks such as precision target tracking (Stevens and Lewis, 1992).

The issue of aircraft control has engendered a novel conundrum for control engineers within the scientific domain. It is widely acknowledged that aircraft dynamics exhibit a profoundly nonlinear nature, manifesting a robust interplay between longitudinal and lateral dynamics. This coupling intensifies markedly as the aircraft undertakes maneuvers at escalated angular rates and heightened angles of attack. Consequently, the development of a dependable controller becomes imperative to counteract such effects, all the while preserving noteworthy resilience against unaccounted-for dynamics and parameter variations. Within the realm of flight control systems, the requisites for optimal performance fluctuate at varying attack angles. Notably, when the angle of attack is low, the primary performance objective lies in attaining impeccable maneuverability. Conversely, at high angles of attack, where the aircraft approaches or enters the stall regime, utmost emphasis must be placed upon preserving flight stability, albeit at the expense of some flight quality compromise. Moreover, in the context of fighter aircraft, the performance requisites may also undergo metamorphosis contingent upon the specifics of flight operations, encompassing factors like speed, attitude, pilot commands, and others (Ijaz et al., 2021). Stability analysis emerges as a momentous phenomenon necessitating comprehensive consideration to achieve the envisaged mission in consonance with the specific aircraft archetype. Noteworthy determinants encompass passenger comfort, the pilot’s command over aircraft manipulation, the meticulous calculation of flight performance, sensor precision, and a myriad of other criteria divulged through an intricate evaluation of the aircraft’s inherent stability. The ever-evolving technological landscape, characterized by amplified technical capacities, heightened maneuverability in temporal domains, and the proliferation of time-varying data such as ammunition and fuel specifics, impels the urgency of performance and efficiency computations alongside rigorous precision assessments. These verifications entail the meticulous utilization of aerodynamic efficiency data stemming from the aircraft’s architectural blueprint, coupled with motion equations predicated upon aerodynamic coefficients, control surface efficacy evaluations, flight performance computations, and a host of other multifaceted considerations (Özcan and Caferov, 2022). Traditionally, flight control systems have been meticulously fashioned through the utilization of mathematical aircraft models, which undergo linearization at multiple operation points, leading to the programming of controller parameters contingent upon prevailing flight conditions (Andrade et al., 2017).

Autopilot systems have demonstrated a significant function in advancing aviation, as they actively enhance navigation protocols, aviation management, and the overall stability and control of the aerial vehicles (Nelson 1998). The inclusion of nonlinear terms in control algorithms introduces intricacy and heightened computational expenses. As a result, PID control algorithms have proven to be effective, owing to their straightforward nature, ease of implementation, and commendable performance across various instances (Kada and Ghazzawi 2011). Due to this nonlinearity, the conventional approach to designing flight control systems involves the utilization of mathematical models of the aircraft that are linearized at different flight conditions. Consequently, the controller parameters or gains are “scheduled” or adjusted based on the flight operating conditions (Vo and Seshagiri, 2008). The F-16 Air Combat Fighter leverages the concept of relaxed static stability (RSS) in the pitch axis, imparting amplified aerodynamic lift and mitigated trim drag. This technological breakthrough constitutes a paramount achievement for the F-16, entailing a state where the aircraft attains equilibrium along the pitch axis during subsonic flight conditions, such that the wing’s center of lift aligns with or precedes the center of gravity. Consequently, the inclusion of a lifting tail becomes imperative. The RSS system exhibits a tendency to swiftly deviate if continuous activation of pitch stability augmentation is not sustained (Ammons, 1978).

In the present work, the progression towards the goals will be methodically executed through a gradual step-by-step approach, ensuring a meticulous achievement of each milestone. Consequently, the study is divided into two distinct sections, with each section comprising two groups. The initial section requires substantial effort and serves as the foundation for subsequent advancements, wherein the tools and knowledge gained from the first section are utilized. In the second section, control systems are designed for the constructed model. Both longitudinal and lateral aspects are individually addressed within each section, resulting in the subdivision of each section into two distinct groups. The dynamics model of system will be constructed by considering its geometrical, mechanical and aerodynamic characteristics of the aircraft. Wherever possible, extracted parameters from literature were listed in tables as much as possible and other needed parameters were calculated employing presented formulations from referred sources. Mathematical model of aircraft was developed in convenient representations to study its dynamic behavior. Then appropriate control systems for both longitudinal and lateral dynamics of modeled system were designed. The PID control method is widely utilized in the aviation industry due to its combination of simplicity and performance. In addition to the classical PID approach, the LQR method, recognized as a modern control technique, has been implemented as an alternative approach within the control systems structure.

2. Flight Dynamics and Modes of Motion; State Space Equations

The study commenced by deriving the flight dynamic equations and identifying the associated parameters. The resulting model will be presented in both state-space and transfer function formats, encompassing both longitudinal and lateral behaviors. The characteristics of these motions and their modes will be elucidated through the examination of the eigenvalues of the state matrix. The state space representation of an aircraft is a mathematical presentation that describes the dynamic behavior of that aircraft using a set of state variables and input-output equations. This representation is commonly used in control engineering and flight dynamics analysis.

The state space model can be articulated in the subsequent manner (Eq 1):

(1)

A          :State Matrix

B          :Input Matrix

C         :Output Matrix

x          :State Vector, Longitudinal; 〖x=[u,w,q,θ]〗^T

x          :State Vector, Lateral; 〖x=[ u,w,q,θ]〗^T

u          :Input Vector

Developing a specific state space representation for an aircraft requires detailed knowledge of the aircraft’s dynamics, aerodynamics, and control systems. The three matrices A to C are typically determined through system identification techniques, simulations, or flight test data analysis. Once the state space model is established, it can be used for various analyses, including stability analysis, control design, and performance evaluation.

The initial aim of this investigatin is to determine the matrices A and B, which are constructed using stability derivatives. The system matrix (state matrix) and input matrix (control matrix) for the longitudinal motion of the aircraft are provided as follows:

(2)

                            (3)

Parameters inside matrices are combination of the stability derivatives. All these elements will be calculated separately.

(4)

               (5)

             (6)

             (7)

                     (8)

                  (9)

                  (10)

                   (11)

                 (12)

                 (13)

                     (14)

                    (15) 

Which the including constants are as follow:

       (16)

(17)

              (18)

              (19)

          (20)

(21)

       (22)

Initial flight condition is supposed to be an altitude of 30,000 ft and flying angle of attack to be five degrees. Starting from atmospheric calculation using given altitude and Mach number. Temperature at an altitude is calculated using equation  and B is 2°C drop in temperature for each 1000 ft ascent. Under ISA conditions, at sea level T₀=15°C so  at desired altitude. Speed of sound calculated using (Eq 23):

(23)

       :specific heat ratio (~1.4 for normal air at S.T.P)

  R        :gas constant (287.26 for air)

  T        :absolute temperature (k)

Then the initial speed of the aircraft can be reached using (Eq 24):

(24)

For calculating air density, following equation (Eq 25) is used:

(25)

Air density at sea level (ρ₀) is 1.225 kg/m³ so at the given altitude ρ=0.46 kg/m³ then dynamic pressure can be obtained using (Eq 26):

(26)

The Table 1 contains the mass and geometric attributes of the aircraft. The mathematical model is based on streamlined high-fidelity data sourced from NASA Langley wind-tunnel experiments performed on a scaled model of the studying airplane (Nguyen et al., 1979). “ ” and “ “will be calculated using model geometry, further explanation is shown in Fig. 1 (Wikimedia drawing). ” ” for biconvex airfoil of tail at M=0.6 achieved 2.86 from literature (Nguyen et al., 1979). It’s assumed that . Calculated data are listed in Table 2. Before starting to calculate stability derivatives, the required nondimensional coefficients can be found in Table 3 from literature (Nguyen et al., 1979) or obtained using given above relations (Eq 16-22).

Table 1. Geometric Characteristics of the Aircraft

Table 2.Geometric Characteristics of the Horizontal Tail

Table 3.Non-Dimensional Derivatives of Longitudinal Stability

Fig. 1. Three view drawing of the F-16 aircraft (Wikimedia drawing), Linear dimensions are in meters

Table 4.Engine Characteristics

Table 5.Longitudinal Stability Derivatives

Now the main parameters, i.e., stability derivatives can be achieved using mentioned relations (Eq 4-15). They are listed in Table 5 and thus forming of state and input matrices is feasible anymore. Two intended matrices (Eq 2&3) are introduced in the following:

(27)

                                   (28)

Previous matrices indicate longitudinal motion of the aircraft. The lateral motion is studied in the succeeding. State matrix “A” and Input matrix “B” are introduced below:

(29)

Parameters inside matrices are a combination of the stability derivatives. All these elements will be calculated separately as previously it’s done for longitudinal section.

        (30)

        (31)

       (32)

      (33)

      (34)

     (35)

     (36)

       (37)

      (38)

   (39)

      (40)

    (41)

   (42)

    (43)

Due to accessibility of all needed non dimensional coefficients from direct data or guessing using available plots in literature (Nguyen et al., 1979). There is no need to recalculate them, these data are listed in Table 6. The vertical tail characteristics are listed in Table 7 to be used for determining some coefficients that have to be calculated.

(44)

         (45)

For simplification of presented relation in (Eq 44) an illustration of the model is shown in Fig.2, an equation can be developed to determine (Eq 44). Using Thales’s theorem for blue triangle in the Fig.2, it’s been resulted:

(46)

By substituting of this term instead of “c” integral term in (Eq 44) will be solved easily:

(47)

Two recently calculated derivatives have been included in Table 6. The elements of the studying matrices have been obtained using  (Eq 30-43), and the resulting parameters are listed in Table 8. The state and input matrices for lateral motion can be formed using relations of (Eq 29). Controllability and observability matrices of “A” are full rank, i.e., state matrix “A” is controllable and observable.

(48)

                               (49)

Table 6.Non-Dimensional Derivatives of Lateral Stability

Fig. 2.Top View of the F-16, to Determine Some Geometric Dimensions (Wikimedia drawing).

Table 7.            Geometric Characteristics of the Vertical Tail

Table 8.Lateral Stability Derivatives

3. Flight Dynamics and Modes of Motion; Transfer Functions

A very useful concept in the analysis and design of control systems is the transfer function. The transfer function gives the relationship between the output of and input to a system. In the case of aircraft dynamics, it specifies the relationship between the motion variables and the control input. To avoid undue mathematical complexity, simpler mathematical models are developed using longitudinal and lateral approximation so that the idea behind various autopilots can be examined (Nelson, 1998). For longitudinal modes, short-period (Eq 50) and phugoid mode (Eq 51) are introduced.

(50)

(51)

Table 9.Constants A, B and C in the (Eq 50) and (Eq 51).

The values of present constants in above equations are listed in Table 9.

Table 10.Transfer Functions of Longitudinal Modes

Substituting model characteristics in contents of Table 9 and using (Eq 50-51), resulted transfer functions are listed in Table 10.

The below relations (Eq 52) also used in (Eq 50) determination.

(52)

The same process is applying for lateral modes; Roll dynamics (Eq 53) and Dutch roll mode (Eq 54).

(53)

(54)

Table 11.Constants A, B and C in the (Eq 54) For Dutch Roll Mode

Presented constants in (Eq 54) detailed in Table 11. Using model characteristics in (Eq 53-54), resulted transfer functions for lateral modes are listed in Table 12.

Table 12.Transfer Functions of Lateral Modes

Eigenvalues of matrix A are roots of the characteristic equation of system and are obtained by solving this equation: . Matrix A for short-period mode is shown in (Eq 55).

 (55)

Eigenvalues for short period mode are resulted “λ₁= -1.1033”, “λ₂= 0.6312”. This mode is including an unstable pole in the right side of the origin. Eigenvalues for phugoid mode are “λ_1,2= -0.0008 ± 0.0767i” which shows almost a stable situation. Matrix A for phugoid mode is shown in (Eq 56).

(56)

Determined eigenvalue for roll mode can be acquired using (Eq 57) that is resulted” ” which is a stable and reasonable root in the left side of the origin.

(57)

For Dutch roll mode we neglected rolling moment equation and stated that depend on sideslip and yawing rate as is shown in (Eq 58). Eigenvalues for Dutch roll mode are “λ_1,2= -0.8272 ± 2.2374i”. These values indicate a stable situation for this mode.

(58)

Furthermore, it is important to mention that the spiral mode, which is considered one of the lateral modes, has not been considered in this study. The corresponding eigenvalue for the spiral mode has been determined to be “ ”. Typically, this value is in close proximity to the origin.

This pole arrangement on s-plane is usual for an aircraft with aft CG (Denieul et al., 2017). As manufacturer says and pilots confirm; this aircraft is very unstable. Moreover, flap effectiveness of this kind of elevator (Stabilator) is one. All these are preconditions for having an agile vehicle. For simulation, engine is supposed to work in half of its maximum thrust ( ) that is about in military thrust range (Nguyen et al., 1979). For lateral states behavior, roll angle is the most affected to aileron input, after that roll rate is sensitive to this control surface. Rudder deflection causes sideslip and yaw moment, since yaw and roll motions inevitably related to each other, as a lateral motion we study, both yaw and roll rates and also sideslip and roll angles are affected enormously

4. Control Modelling; Autopilot Design

In previous section it was discussed about two types of aircraft stability and achieved related mathematical model in both state space and transfer function representations. In this section the objective is embedding an appropriate controller in these models and reaching desired and reasonable states. A simple pitch control autopilot has been developed to maintain longitudinal stability of the model; using conventional PID and LQR, trying to hold pitch angle close to desired input as much as possible without any noise in final signal. By using inner rate feedback loop we are able to stabilize system. The model’s block diagram is shown in Fig. 3.

In the event that the PID controller parameters (namely, the gains of the proportional, integral, and derivative terms) are inappropriately selected, the controlled process input may become unstable, exhibiting a diverging output, either with or without oscillations, and limited solely by saturation or mechanical failure. The underlying cause of instability lies in an excessive gain, especially in the presence of significant lag. Generally, response stabilization is imperative, ensuring the absence of oscillations under all combinations of process conditions and set points, although, on certain occasions, marginal stability (bounded oscillation) may be deemed acceptable or desirable. Introducing a PID controller to an F-16 aircraft will lead to an enhanced steady-state solution, reduced rise time, and an improved transient response, thus yielding greater accuracy in its performance (Sayegh, 2014).

Fig. 3. Pitch/Roll autopilot block diagram in state space form using PID controller

It can be demonstrated that if the system is state-controllable, it becomes feasible to devise a linear control law that can achieve any desired closed-loop eigenvalue structure. For a single-input system, the control law is represented by (Eq 59), where “η” denotes the control input, and “k” is a column matrix or vector of unknown gains. The Bass-Gura method offers a straightforward approach to determine the gains required for a specific eigenvalue structure. The plant matrix, in general, may not be in the companion form. If the system is not in the companion form, we can employ a transformation using (Eq 60), wherein “V” denotes the controllability test matrix, “W” forms a triangular matrix based on the coefficients of the open-loop characteristic equation, and “ “ and “a” represent the coefficients of the desired closed-loop characteristic equation and the coefficients of the open-loop plant matrix characteristic equation, respectively (Nelson, 1998). This parameters for studying model are listed in Table  13. By incorporating the “k” parameter into the block diagram illustrated in Fig. 3, and leveraging the inner loop(s) and PID controller, a robust methodology can be established to address instability and achieve desirable outcomes within the control systems architecture.

                    (59)

(60)

Table 13.Performance values of PID

Fig. 4. Pitch/Roll autopilot block diagram using LQR controller

LQR controller as a modern control approach is another strategy which is used here to build the pitch autopilot. This autopilot’s block diagram including state-space matrices and controller gains is shown in Fig. 4.

Results from two approaches are shown in Fig. 5. For longitudinal pitch control. Investigating the results for Pitch revealed that rise time from PID method is about 0.1 second, settling time for two percent of final value is around 1 second, with 14 percent overshoot. Without accounting on inner loop, the system at most could be marginally stable with lots of disturbing oscillations and here the importance of the inner loop in increasing damping of the system is obvious. From LQR results rise time was around 2.3 seconds and settling time around 4 seconds with no overshoot, which is too slow.

Fig. 5. Pitch control using top) PID down) LQR controller

In the F16 aircraft, both the aileron and the stabilator, which functions as a differential tail, contribute to roll and lateral motion to enhance maneuverability. When the aileron is deflected by 4 degrees, the differential tail deflects approximately 1 degree. However, the involvement of the differential tail in the system has not been considered in this study. The roll attitude autopilot design is presented in Fig. 3, while the control block diagram for the roll autopilot, utilizing the LQR approach, is depicted in Fig. 4. The system’s response to these configurations is illustrated in Fig. 6.

Fig. 6. Roll control using top) PID down) LQR controller

5. Results and Discussion

PID control is a classical control method widely used in various engineering applications, including aircraft control systems. It operates based on error feedback, continuously adjusting the control inputs (control surfaces) in proportion to the contrast between the desired and current states. On the other hand LQR is a modern control technique based on optimal control theory. It considers a cost function that quantifies the system’s performance and aims to minimize it. LQR is capable of providing an optimal control law that balances the trade-off between control effort and performance. The choice between PID and LQR control for the F-16 pitch and roll autopilot design depends on several factors, including performance requirements, implementation complexity, and available resources. LQR generally provides superior performance compared to PID control in terms of precision, stability, and ability to handle disturbances. LQR’s optimal control formulation allows for fine-tuning control efforts to achieve desired performance. PID control is relatively straightforward to implement and tune. It has been used successfully in many control systems, including older aircraft models. However, it may require more tuning effort to achieve desired performance. LQR is more adaptable to various operating conditions and can handle a broader range of aircraft dynamics. As aircraft become more advanced and incorporate additional complexities, LQR’s flexibility becomes advantageous. PID control is computationally lighter and may be preferred for systems with limited processing power. LQR, being an optimal control approach, might require more computational resources. Finally, the implementation complexity and available resources should also be taken into consideration when making the final design decision.

Whitin this study, the PID controller yields a rise time of approximately 0.2 seconds and a settling time of about 3 seconds for two percent of the final value, accompanied by an overshoot of 8 percent. On the other hand, employing the LQR approach results in a rise time of around 2 seconds and a settling time of approximately 4 seconds, with no overshoot. The second approach demonstrates the absence of overshoot; however, it is slower in comparison. In contrast, the first approach exhibits a slight overshoot, which is inevitable, but offers significantly faster performance. The detailed results obtained from the simulations are presented in Table 14. The outcomes for both the pitch and roll autopilots are deemed satisfactory and reasonable when compared to the existing literature (Ahmed et al., 2019), (Stachowiak and Bosworth, 2004) and (Lu and Wu, 2005). Here, considerably superior outcomes are achieved in comparison to those presented in the initial literature. Moreover, it is important to acknowledge that a variety of results may be observed within an acceptable range, attributable to the collection of information from diverse sources under distinct operating conditions. Furthermore, by altering the Q and R matrices in the LQR strategy, superior results can be attained. Nevertheless, it is preferable to adhere to the main approach to avoid the complexities of model instability and ensure the validity of the simulation data.

6. Conclusions

Within this study, state-space representations were constructed for both longitudinal and lateral dynamics of the system. Additionally, transfer functions were derived to describe the dynamic modes. Subsequently, two autopilots were designed for pitch and roll control, employing both PID and LQR methods. The PID technique has demonstrated its remarkable effectiveness in the design of autopilots, surpassing the performance achieved by the LQR control system. This superiority is evident in Table 14, where the results of the PID approach outperform those of the LQR control system across various performance criteria under this typical simulation conditions. 

Table 14.Performance values of PID and LQR

In future research endeavors, it is recommended to enhance the realism of the study by incorporating the dynamics of actuators through the utilization of appropriate transfer functions. Additionally, the dynamic representations obtained can be employed to design additional autopilot systems using both the PID and LQR approaches. While it can be acknowledged that modern control theories, such as “Adaptive Control”, “Model Predictive Control (MPC)”, “Sliding Mode Control (SMC)”, “H-infinity Control”, “Fuzzy Logic Control” and “Neural Networks” are known to exhibit robust performance in autopilot design and contribute to the overall effectiveness and safety of aircraft operations, The primary aim of this study was to establish a comprehensive understanding of the principles underlying analytical modelling of aircraft dynamics and the analysis of its control systems. The main focus was to take the initial step in this field of study with precision, thereby establishing a solid groundwork for future research and exploration.

Nomenclature

S.T.P      : Standard Temperature and Pressure

b             : Wing Span

S             : Wing Area

c             : Mean Aerodynamic Chord 

CG          : Centre of Gravity

AR          : Wing Aspect Ratio

HT         : Horizontal Tail

rc           : Root Chord

tc           : Tip Chord

b_t            : Tail Span

l_t             : Distance Between CG and Aerodynamic Centre of HT

S_t            : Horizontal Tail Area

V_H           : Tail Volume Ratio

d𝜀/dα   : Downwash Change

𝜂             : Flap Effectiveness

S_v           : Vertical Tail Area

S_R           : Rudder Area

τ_r            : Rudder Effectiveness

τ_a            : Flaperon Effectiveness

S_f            : Flaperon Area

K            : Empirical Factor

VT          : Vertical Tail

l_v             : Horizontal Range Between Aerodynamic Centre of               VT and CG

Z_v           :Vertical Range Between Aerodynamic Centre of VT               and CG

CRediT Author Statement

Masoud Norouzi: Methodology, Software, Data Curation, Validation, Formal Analysis, Writing-Original Draft, Visualization, Writing-Review & Editing.
Elbrus M. Caferov: Conceptualization, Project Administration, Supervision, Resources.

References

Ahmed, W., Li, Z., Maqsood, H. and Anwar, B., (2019). System Modelling and Controller Design for Lateral and Longitudinal Motion of F-16. Automation, Control and Intelligent Systems. Vol. 7, No. 1, 2019, pp. 39-45. doi: 10.11648/j.acis.20190701.15

Ammons, E.E., (1978). F-16 Flight Control System Redundancy Concepts, General Dynamics, Fort Worth Division, Fort Worth, Texas, USA

Andrade, J.P.P, Campos, V.A.F., Potts, A.S. and Garcia, C., (2017). Damping Improvement of a F-16 Aircraft through Linear Matrix Inequalities. International Federation of Automatic Control (IFAC) Hosting by Elsevier Ltd., IFAC PapersOnLine 50-1 (2017) 3947–3952

Denieul, Y., Guibé, J.B., Alazard, D., Toussaint, C., Taquin, G., (2017), Multicontrol Surface Optimization for Blended Wing-Body Under Handling Quality Constraints. Journal of Aircraft, American Institute of Aeronautics and Astronautics ,pp.1-14. DOI: 10.2514/1.C034268

Ijaz, S., Fuyang, C., Hamayun, M.T. and Anwaar, H. (2021). Adaptive integral-sliding-mode control strategy for maneuvering control of F16 aircraft subject to aerodynamic uncertainty. Applied Mathematics and Computation 402 (2021) 126053, www.elsevier.com/locate/amc

Kada, B. and Ghazzawi, Y., (2011), Robust PID Controller Design for an UAV Flight Control System. World Congress on Engineering and Computer Science (WCECS), San Francisco, USA, ISBN: 978-988-19251-7-6

Li, B. and Wu, F., (2005). Probabilistic Robust Control Design for An F-16 Aircraft. AIAA Guidance, Navigation, and Control Conference and Exhibit, San Francisco, California

Nelson, R.C., (1998). Flight Stability and Automatic Control, WCB/McGraw-Hill, ISBN 0-07-046273-9

Nguyen, L.T., Ogburn, M.E, Gilbert, W.P., Kibler, K.S., Brown, P.W. and Deal, P.L., (1979). Simulator Study of Stall/Post-Stall Characteristics of a Fighter Airplane With Relaxed Longitudinal Static Stability, NASA Technical Paper 1538, Langley Research Center, Hampton, Virginia, available at: https://www.cs.odu.edu/~mln/ltrs-pdfs/NASA-79-tp1538.pdf (accessed 20 April 2022)

Özcan, A.B. and Caferov, E., (2022). Frequency Domain Analysis of F-16 Aircraft in a Variety of Flight Conditions. International Journal of Aviation Science and Technology, Volume 3, Issue 1, (2022), 21-34, DOI: 10.23890/IJAST.vm03is01.0103

Reichert, G., (1992/1993). Flugmechanik III: Flugeigenschaftskriterien, Elastisches Flugzeug und Aktive Steuerung, Vorlesungsmanuskript, Institut fuer Flugmechanik des TU Braunschweig, Germany

Sayegh, Z.E. and Deghidy, A., (2014). Auto Pilot Design for F-16. Technical Report, A Project Submitted to the Graduate Faculty of The University of Concordia, Montreal, Quebec, Canada, DOI:10.13140/RG.2.2.36709.91362,  https://www.researchgate.net/publication/325450374

Stachowiak, S.J. and Bosworth, J.T, (2004). Flight Test Results for the F-16XL With a Digital Flight Control System. NASA/TP-2004-212046, NASA Dryden Flight Research Center Edwards, California, USA

Stevens, B.L. and Lewis, F.L., (1992). Aircraft Control and Simulation, A Wiley-Interscience Publication, John Wiley & Sons, Inc.

Vo, H. and Seshagiri, S., (2008). Robust Control of F-16 Lateral Dynamics. International Journal of Mechanical, Industrial and Aerospace Engineering, DOI: 10.1109/IECON.2008.4757977, Source: IEEE Xplore

Wikimedia drawing, A 3-view line drawing of the General Dynamics F-16 Fighting Falcon, Page URL: https://commons.wikimedia.org/wiki/File:General_Dynamics_F-16_Fighting_Falcon_3-view_line_drawing.svg