Review on System Identification and Mathematical Modeling of Flapping Wing Micro-Aerial Vehicles

Featured Application

The paper provides a comprehensive review on various system identification and modeling techniques for flapping wing micro-aerial vehicles (FWMAVs) that will be useful for researchers in order to obtain mathematical modeling of FWMAVs.

Abstract

This paper presents a thorough review on the system identification techniques applied to flapping wing micro air vehicles (FWMAVs). The main advantage of this work is to provide a solid background and domain knowledge of system identification for further investigations in the field of FWMAVs. In the system identification context, the flapping wing systems are first categorized into tailed and tailless MAVs. The most recent developments related to such systems are reported. The system identification techniques used for FWMAVs can be classified into time-response based identification, frequency-response based identification, and the computational fluid-dynamics based computation. In the system identification scenario, least mean square estimation is used for a beetle mimicking system recognition. In the end, this review work is concluded and some recommendations for the researchers working in this area are presented.

1. Introduction

The fascination toward the flight of natural fliers has inspired the development of flapping-wing micro air vehicles (FWMAVs). Natural fliers, particularly of insect types, offer interesting flight mechanisms and aerodynamics that have not been completely understood until today. The currently estimated more than ten million insect’s species [1] include flies, which change their flight conditions based on some specific arrangements of their wing strokes [2,3], offer different adjustments of the wings that make the insects capable of flying with extreme maneuverability and impressive stability [1,4,5] compared to the conventional fixed wing aircrafts [6,7]. Some characteristics like their dexterity, broad flight envelope, robust behaviors in different situations at very low Reynolds numbers, hovering ability, the ability of quick transition to forward flight, abrupt maneuvers, make FWMAVs potentially beneficial in this era. Insect fliers generally use elastic elements in their flight muscles as energy storage and thus, minimize their inertial power, thus enhancing their power loading [8]. This results in a more efficient flight, which is another advantage of flapping-wing flight over rotary wings [9]. Exhaustive applications of such MAVs can be found, for example, as a reconnaissance platform for search and rescue as well as military operations. Fitted with a camera, FWMAVs can be deployed in risky and hazardous conditions such as after natural catastrophes to quickly determine and monitor the situation, for finding routes, identify potential dangers, and patrol border lines for security purposes as well as for civil applications such as the remote inspection of pipes and power lines, inspection of crop fields, and the conditions for fertilization and harvesting [10,11,12,13].

Based on the current trends, it is envisioned that FWMAVs will dwell in the gaps between the rotary and fixed wing aircraft because of their aforementioned appealing characteristics. In this context, micro- and nano-flapping wing fliers [10,11,12,13,14,15,16,17,18,19,20] are also believed to provide a wide number of advantages with considerably low cost and flight efficiency. As shown in [21], insect flights show higher power loading than multi-rotor drones. However, for FWMAVs, the small size of the sensors and mass restricts the onboard processing. Furthermore, depending on the overall size of the MAVs, the limited choice of small size onboard power supplies often means that the flight duration is a relatively short period in comparison with other types of aircraft. Thus, at the present stage of development, such efficiency has yet to be achieved. Therefore, energy efficient flight, attitude determination, and flight control of FWMAVs are still the main focus of researchers, scientists, and engineers [13,16].

FWMAVs are primarily recognized and developed into two classes (i.e., tailless and tailed). The tailless designs usually mimic the insects, which, consequently, provide high dynamic maneuvers. These configurations, in general, require active stabilization because they can be found passively unstable in forward flight and hovering [22]. In such developments, the nano-hummingbird [7] and robotic bee, known as RoboBee [10], have shown actual flight stability. On the other hand, the tailed design is found to be similar to configurations in birds, which are known for being passively stable, and therefore do not require active stabilization [15,23]. This kind of configuration provides the ease to develop and implement higher-level problems like altitude stabilization [24,25] and obstacle avoidance [25,26,27]. Figure 1 displays some of the well-known flapping-wing micro air vehicles that have been developed to date.

It is witnessed in Figure 1 that these FWMAVs differ in design as well as in its flying mechanism. The FWMAVs (Figure 1a,b,f,h) follow very similar flight mechanics as birds. This class utilizes the surface control mechanism for different maneuvers like pitching, rolling, yawing, and hovering. The surface motion control resembles the fixed wing aircraft strategies [7], whereas the FWMAVs in (Figure 1c,d,e,g) mimic the insects flying mechanism to experience different flight maneuvers [11,14,15,16,19,20,28,29]. The autonomous flight of the first category is relatively easier compared to that of the second class of FWMAVs. A clear understanding of insect flight depends on subtle details that might be easily overlooked in otherwise thorough theoretical or experimental analyses. Recently, researchers have benefited greatly from the use of high-speed cameras for capturing wing kinematics, new strategies such as digital particle image velocimetry (DPIV) to quantify flows, and powerful computers for further computational analysis. With the help of these methodologies, one can develop more accurate and rigorous models of these insects while having a fewer number of assumptions in developing wing kinematics, aerodynamics, and flight dynamics.

Despite the challenges faced in conducting experiments to quantify the wing motions of several flying insects due to their small size and high wing beat frequencies, insect fliers’ flapping kinematics, particularly that of free hovering flight and free forward flight, have been well-researched for more than 50 years. Early attempts to capture free-flight wing kinematics such as Ellington’s comprehensive and influential survey [30] relied primarily on a single-image high-speed cine. Although quite informative, especially because digital imaging technology continues to offer exceptional spatial resolution, single-view techniques cannot provide an accurate time course of the angle of attack of the wings. More recent methods have employed high-speed videography [31], which offers greater light sensitivity and ease of use, albeit at the cost of image resolution. A further consideration is that insects rely extensively on visual feedback, and hence, care must be taken to ensure that lighting conditions do not significantly impair an insect’s behavior. Even more challenging than capturing wing motions in 3-D is measuring the time course of aerodynamic forces during the stroke. At best, flight forces have been measured on the body of the insect rather than its wings, making it very difficult to separate the inertial forces from the aerodynamic forces generated by each wing [32,33,34,35,36]. In addition, unlike in free flight conditions, tethering can alter the wing motion, and thus, modify the generated forces. Researchers have overcome these limitations using two strategies. The first method involves constructing dynamically scaled models on which it is easier to directly measure aerodynamic forces and visualize flows [37,38,39,40,41,42,43]. A second approach is to construct computational fluid dynamic simulations of flapping wings [44,45,46,47,48]. The power of both approaches, however, depends critically on accurate knowledge of wing kinematics.

In general, to understand the mechanics of flapping wing flight, three basic aspects need to be well understood: flapping kinematics, aerodynamics, and flight dynamics (a detailed summary can be found in [49,50,51]). Several attempts have been made to model FWMAVs based on the knowledge of those three aspects of natural fliers such as those in [52], which focused on the first principle approach of modeling that includes flight stability and control. However, since FWMAVs depend on the flapping mechanism, the actuation system for generating the aerodynamic forces and moments as well as the sensors and control dynamics, mathematical models of FWMAVs solely based on the first principles may not accurately describe the actual dynamics of the MAVs. The review provided in this paper is intended to give insights into system identification approaches as an alternative way in obtaining mathematical models of various FWMAVs based on experimental data. Such approaches differ from the conventional modeling approaches through physics (i.e., modeling through the first principles). With system identification approaches, the richness of experimental or flight data can be exploited to obtain the flight dynamics and aerodynamics characteristics of FWMAVs. The models obtained through such approaches can be used to validate mathematical models (including those obtained through the first principles) and predict the effect of any design changes.

System identification uses mathematical modeling techniques, which are based on observed input–output data. This process was first coined by Lofti Zadeh in the late 1950s [53]. It has been used successfully to identify models of various complex systems such as the cardiovascular system [54], aircraft [55,56], and underwater vehicles [57]. This type of modeling approach is beneficial in developing realistic mathematical models of complex systems based on experimental data for prediction, simulation, system design, and/or control [58]. System identification is usually conducted via (1) the Black box approach in which system identification is carried out by fitting the system’s input–output response with some flight data, or (2) the Grey box approach (see [59] for more details), which is carried out by first developing the structure of models mathematically, then the required parameters in the models are identified. In both approaches, the overall models are validated against the experimental data. The typical steps involved in system identification are shown in Figure 2.

In this paper, the most prominent strategies relevant to system identification of the flight dynamics and aerodynamics of FWMAVs as well as natural fliers were overviewed. This paper is organized as follows. A general description of FWMAVs is presented and some major works done so far are highlighted in the Introduction. In Section 2, the background of the flight of FWMAVs is presented comprehensively and recent developments in this area are highlighted. Section 3 contains a brief discussion on system identification approaches, which then leads to the main review that presents different strategies in system ID adapted for modeling the FWMAVs’ dynamics. The article is presents our concluding remarks and some recommendations in Section 5, where the authors present their own proposed strategy for an online system’s stability derivative estimation that is potentially advantageous in obtaining the flight dynamics of any existing FWMAVs as future works.

2. Micro-Aerial Vehicle Flight

Current developments in flapping wing MAVs, which are inspired by natural fliers, can be categorized into two classes: bird (or vertebrate) flight and insect flight [60]. Although in general, both types of natural fliers observe three types of flight, which are named as gliding flight (mostly found in vertebrate flight), forward flight, and hovering [12,61,62,63,64,65,66], there are differences in their aerodynamics (see Figure 3). In bird flight, the wings generally have two degrees of freedom (i.e., the flapping and symmetric twisting of the wings downward and upward, i.e., pronation and supination occur respectively) to change the angle of attack and the lift. However, in the second class, the wings follow three degrees of motion (i.e., the main flapping motion is defined along the body axis with a certain stroke plane, while flapping the wing exhibits pitching motion along its root and experiences deviation from the stroke plane). In other words, while flapping, the wings follow an eight like pattern, which is, by definition, a two-degree motion [52]. Consequently, a three degree of motion occurs. The wing strokes are divided into upstroke and downstroke that are further categorized into normal rotation, advanced rotation, and delayed rotation. The wing rotation is termed as normal if the rotation occurs at the end of each half stroke, the rotation is advanced if it occurs before the end of the half stroke and after the mid stroke, and delayed rotation occurs after the half stroke and before the mid stroke. Furthermore, the normal rotation is, once again, recognized as water treading and normal hovering (for more details see [60,67,68]). Detailed comparison between the vertebrate and insect flights can be found in [69,70].

In the context of understanding the aerodynamics and wing kinematics, the lift and power requirement analysis related to hovering natural fliers (e.g., hummingbird and Drosophila family) was conducted in seminal works [73,74], which was further refined in the work of Ellington [30,75,76,77,78,79]. These pioneering works were further developed for other species, which include the design of some morphological parameters and wing kinematics (e.g., hawkmoths [80,81], dragonflies [82,83], bumblebee [84,85], robotic bee [10,20], and the beetle [11]. The flight characteristics are modeled in terms of aerodynamics and wing kinematic models. In the existing literature, two main strategies have been extensively studied for the development of the aerodynamic models (i.e., blade element methodology (BEM) and Navier–Stokes equation) through the computational fluid dynamics (CFD) strategy. The blade element strategy based aerodynamic models was developed in [61,68], where they incorporated the unsteady and added mass effects. The progress in the aerodynamics modeling in the last decades is comprehensively presented in [52,86,87,88,89,90,91,92].

In recent developments, [93,94] tried to abridge the theoretical and experimental strategies and to develop a standard linear aerodynamic model for DelFly II. The developments were mainly based on the standard aircraft system identification techniques and their extension to FWMAVs, which can be further used for simulation and autonomous full flight envelope. In [95,96], a two-dimensional flapping wing was considered, and the aerodynamic performance in asymmetric stroke in hovering and forward flight was numerically studied. In addition, the analysis of the effects of the asymmetry of the stroke on the aerodynamic forces and flow structures was carried out. The authors in [97] carried out a BEM for the development of the fluid–structure interaction model by synthesizing it with a structural model. The benefits of this strategy were that it not only improved the reliability of the BEM, but it also enabled various aerodynamic studies (e.g., [98,99,100,101,102]). Note that the BEM method is widely utilized due to its simplicity and robustness, however, its accuracy has been questioned. In contrast, computational simulations based on the solution of Navier–Stokes equations [103,104,105] have been utilized because they provide, with relatively high accuracy and precision, the detailed flow structures, time-dependent aerodynamic forces, and power requirements. However, this strategy needs high-fidelity modeling and a long experiment processing time, which in turn diminishes its usage in practice.

A CFD-informed quasi-steady model based on a hybrid of high- and low-fidelity models to reduce the aforementioned complexity was presented in [106]. The systematic development regarding the beetle’s and DelFly II’s aerodynamics can be found in [107,108,109], respectively. In order to update the research carried out in the area of FWMAVs, a comprehensive bibliometric review can be found in [110]. The flight stability of different species of natural fliers, particularly the longitudinal stability, was studied successfully by a number of researchers, for example, the stability of desert locust Schistocerca Gregaria was thoroughly presented in [86], the nano-humming bird was focused on in [12], insects mimicking the FWMAV’s initial stability was addressed in [111], and bumblebee stability was presented in [112]. In general, the eigenvalue-based analysis resulted in three modes common to all fliers, namely, the unstable oscillatory mode, stable fast subsidence mode, and stable slow subsidence mode. Almost all mathematical models developed in the aforementioned studies fall in the category of modeling by the first principles. System identification approaches have hardly been used. Therefore, the forthcoming section will present, in detail, the techniques devised for system identification approaches, particularly, in obtaining system stability derivatives and control-oriented dynamic models of FWMAVs and natural flapping fliers.

3. System Identification and Mathematical Modeling of Flapping Wing Micro Air Vehicles (FWMAVs)

The ultimate objective, in the development of FWMAVs, is securing an autonomous flight in order to utilize them for the different tasks mentioned earlier. Therefore, a first step toward the establishment of an autonomous flight is the modeling of the flapping kinematics, aerodynamics, and flight dynamics of FWMAVs. This task is particularly challenging due to inherent instability, nonlinear, time-varying, and periodic nature of the vehicles’ dynamics, especially during dynamic aerial maneuvers such as sharp and abrupt turns, obstacle avoidance, and mid-air prey capturing. The aerodynamics and flight dynamics of such maneuvers, even for natural fliers, have not been fully uncovered yet. Additionally, on-board sensors, flapping mechanism, and actuation systems may add complexities to the overall dynamics of the MAVs. Furthermore, for insect-mimicking FWMAVs, due to their small size, the selection of hardware for their avionics and flight control systems are limited. However, in this review, we focused on system identification approaches in order to model the FWMAVs without actually going through the details on hardware selection.

In general, two complementary strategies are applied to model FWMAVs. These strategies are named as bottom-up and top-down [94]. The bottom-up approach mainly focuses on the fundamental approaches to explain and model the behavior and interaction of air motion with the flapping wings. This is essentially modeling with the first principles based on the well-studied aerodynamics and wing kinematics of natural flapping fliers (see [50] for a detailed summary). The practical system development of bio-inspired FWMAVs remains the later part in such developments. On the other hand, in the top-down approach, the construction and developments of the real MAVs take place prior to the theoretical developments. In order to obtain the mathematical model(s) of the FWMAVs, the first step is the system identification, which means a strategy for acquiring a known mathematical model (for instance, aerodynamic modeling, system stability derivatives, and system’s control derivatives) of the system while utilizing the available measured data of the system’s input and output signals produced by the MAVs. These available signals may be measured either in time domain or in frequency domain. The main advantages of performing system identification in modeling are that the process uses real-time data and can take into account design changes. The latter is particularly useful for control system design. Through the process, the value of specific design parameters and their effects on the dynamics of MAVs can be determined, and the structure of the mathematical models obtained is not unique, which gives the freedom to the engineers/researchers to select the models that suit the objectives of the modeling. For examples, these models may be suitable for (hardware-in-the-loop) simulations, analysis, and/or control.

With regard to this stage, the FWMAVs and their flight have been enlightened in a general way, and now the objective is to present the system identification in detail. It is generally conducted either for equivalent linear systems or core nonlinear systems. This may further be classified into time-domain and frequency-domain based identification for tailed and tailless FWMAVs. It is obvious that the dynamic modeling of FWMAVs always considers the mechanism kinematics of linkages, which connect the motors and actuators to FWs [113], aerodynamic forces [114], and nonlinear phenomena like saturation and hysteresis [115]. Consequently, it results in different equations and realizations of the dynamics.

The linear modeling and identification of fixed wing aircraft were extensively studied in [55,114]. The identification strategies for modeling FWMAVs or natural fliers vary from system to system because of their varying kinematics and flight dynamics. The review paper by [52] highlighted that the inertia effects of the wings on the main body were invoked in very few models of the flapping wing systems. The majority of the flapping angles of FWs remains greater than 20 degrees (see, for instance, [116]). Therefore, the small angle assumption (i.e., $\sin \left(\theta \right)\approx \theta$) is no longer valid, and consequently, the linearized dynamic model becomes inaccurate. The dynamic models of the FWs are generally highly nonlinear and time varying due to the aerodynamic forces on the wings that are typically predicted via the lift and drag coefficients of the quasi-steady aerodynamic model. These coefficients depend on the angle of attack, which may vary significantly during a single stroke. Various nonlinear modeling approaches have focused on the addition of a spring in the redesign and optimization of the elastic elements like actuators, wings, and linkages to further improve the capability of the system to be driven at resonances [16,17,18,19]. Another strategy [20] considered wing, actuator, and transmission components, and developed an integrated nonlinear model, which considers all the changes made in the wing shape and linkage system. The change in the wing shape significantly affects the inertial properties and the aerodynamic damping forces, whereas the variation in the linkage mechanism affects the control input channels. Some useful linear models have also been developed that account for the nonlinear geometric and aerodynamic components [117,118]. The developed nonlinear model sufficiently predicts the behavior of the FWMAVs with a considerable operating range. However, the accuracy in the model reduces when operated at high frequencies.

The FWMAVs, as mentioned earlier, mimic either insect or bird flight. They are stabilized via active and passive stabilization strategies in different flight envelopes [10,12,30,31,32] through manually tuned controllers such as proportional-derivative (PD) and/or proportional-integral-derivative (PID) controllers. However, they are observed to be sensitive to external disturbances. Thus, it becomes quite challenging to keep them stable in broad flight envelopes. One reason for such difficulties is the absence of a sufficiently accurate nonlinear model that can describe the unsteady aerodynamics of the FWMAVs’ flapping motions [119,120] and the added inertia effect contribution to the MAVs dynamics [121]. In this context, a linear time invariant (LTI) model was developed in [121,122,123,124] via the use of standard aircraft equations. Another approach, known as the multi-body approach, was introduced. However, such an approach for a time varying nonlinear model [121,125] is not suitable for onboard processing. All the aforementioned approaches have significant shortcomings such as a lack of practical validation, limited only to computational environment and the repeatability of the data. Therefore, it is demanding to overcome such modeling difficulties in order to fully explore the capabilities of these fliers.

In the subsequent subsections, we will focus on system identification strategies that have been applied for identifying parameters and mathematical models of FWMAVs. In presenting the review, we classify natural fliers and FWMAVs into two categories, the tailed and tailless fliers/FWMAVs. Various system identification techniques used in identifying the model parameters of fliers/FWMAVs in these two categories will be discussed.

3.1. Tailed Flapping Wing Micro-Air Vehicles

In this section, we intend to present the details of system identification techniques for tailed micro air vehicles (i.e., bird like micro air vehicles). In [126,127], a flight test of a 0.45 kg ornithopter was flown in a straight and level mean flight. A visual tracking strategy was employed to provide measurements of the states, which benefited with low-noise, minimally-invasive measurements. The spatial distributions of kinematics over a wing stroke were investigated while using a multi-body model of the flight dynamics. The models of lift, thrust, and the coefficients of pitching moment were extracted via stepwise regression accompanied by least mean square estimation. Both the frequency and time domain system identification strategies were invoked, which resulted in quite close estimates. The identified models can be used for feedback control purposes. A grey box approach was proposed for the system identification of the DelFly II in [93,94,128], which by definition belongs to the top-down approach. The purpose of their work was to bridge the existing gap between the theoretical and experimental strategies. Another core objective was the development of benchmark aerodynamic models for the DelFly II based on aircraft system identification techniques. In addition, the extension of these techniques to the flapping wing system for autonomous flight was the key objective. The identification was mainly done via the least mean square estimation and extended Kalman filter (EKF), which are highlighted in the subsequent study.

3.1.1. System Identification of DelFly II Based on Least Mean Square Estimation

The Delfly II, shown in Figure 4, is a bioinspired ornithopter, which is configured with four flapping wings and an inverted ‘T’ tail. Its weight is 17 g and is capable of hovering flight as well as forward flight at the speed of 8$m{s}^{-1}$ [93]. The experimental flight was conducted in an experimental chamber, which was equipped with a tracking system capable of recording the position of a small reflective marker at 200 Hz frequency. The position and altitude were tracked via the use of eight markers placed on its structure. A set of longitudinal elevator input with step, doublets, and triplets with a reference time of 0.3 s was programmed into its autopilot to acquire a full dynamic response over a set of frequencies. The framework for the system identification was defined as follows. The first step toward the system identification was observation/estimation of the state variables of the system from the provided position data. This process was named as flight path reconstruction in [129]. Then, these estimated and available data were used to reconstruct the indirectly measurable states such as angle of attack and refinement of the estimated states. The whole reconstruction was carried out in two reference frames. The tracking system frame was named as the inertial reference frame and was indicated with a subscript I at the respective coordinate axes so that the positive ${\overrightarrow{Z}}_I$ was pointing upward. The frame fixed with body was identified by a subscript b at the respective frame. The positive ${\overrightarrow{x}}_b$ was described along the forward path, ${\overrightarrow{z}}_b=normalized ({\overrightarrow{i}}_b\times {\overrightarrow{W}}_{LR})$, where ${\overrightarrow{i}}_b$ is the unit vector along ${\overrightarrow{x}}_b$ and ${\overrightarrow{W}}_{LR}$ is a unit vector from the left to right wing marker vector. Furthermore, ${\overrightarrow{y}}_b$ was defined along the $({\overrightarrow{x}}_b\times {\overrightarrow{z}}_b$). Consequently, two reference frames were defined for the analysis and reconstruction of the states. The Euler angles were considered according to the standard literature approach of yaw (ψ), pitch (θ), and roll (ϕ). Any vector ${\overrightarrow{v}}_I$ in the inertial reference frame can be defined in the body frame coordinated via the following relation.

$${\overrightarrow{v}}_b=R_{bI}{\overrightarrow{v}}_I,$$

(1)

where $R_{bI}$ is the standard rotation matrix [6,7].

It is witnessed that in the near hover configuration of the DelFly, the yaw angle (ψ) and roll angle (ϕ) observed an abrupt variation of (${180}^{°}$). The Euler angles’ yaw (ψ), pitch (θ), and roll (ϕ) were converted into new coordinates (${\varphi }_r$), (${\theta }_p$), and (${\psi }_y$), which can be above (${90}^{°}$) and to keep them valid for the aircraft equation of motion [6,7]. In addition, the dynamic equations were allowed, via the above coordinated conversion, to remain valid for the pitch angles in the interval (−π/2, 3π/2). From the available data, state reconstruction including the angular rates, angular accelerations, linear velocities, and linear acceleration was accomplished via the discrete time approach, and the effects of high oscillatory modes were reduced by using a 3rd order Butterworth low pass filter with the cut-off frequency of 30 Hz. The angular velocities were recovered via the following equation:

$${\left[\begin{array}{c}p\\ q\\ r\end{array}\right]}_t=\left[\begin{array}{c}{\varphi }_{rt}-{\varphi }_{rt-1}\\ {\theta }_{pt}-{\theta }_{pt-1}\\ {\psi }_{yt}-{\psi }_{yt-1}\end{array}\right]\times \tau =\left[\begin{array}{c}\mathrm{∆}{\varphi }_r\\ \mathrm{∆}{\theta }_p\\ \mathrm{∆}{\psi }_y\end{array}\right]\times \tau .$$

(2)

To obtain the full nonlinear dynamic equations for the entire flight envelope of the aircraft is quite challenging. In this work, the nonlinear model was approximated with a linear model by following [55,130,131] for sufficient small excitations in the neighborhood of the stationary trimmed regime. It was assumed that the dynamics of the flyer could be modeled as Newton’s Euler equations of motion while further assuming that the aircraft body is rigid with constant mass and no changing inertia, subject to no flapping in the absence of wind. Note that the flapping was modeled as a thrust force. The aerodynamic model structure was chosen as follows:

$$F_i=C_{F_0}+{\displaystyle \sum }_{S=1}^n{C}_{F_s}.S; S=\left\{{\varphi }_r, {\theta }_p,{\psi }_y, u, v, w, p, q,r, \alpha ,\beta ,\gamma ,{\delta }_f, {\delta }_e, {\delta }_r\right\}$$

(3)

where $C_{F_0}$ is the affine coefficient; $C_{F_s}$ the states’ coefficient vector; and S is the state vector. These coefficients as well as the states were estimated via the discrete time differentiation. Equation (3) represents the full model structure so that the aerodynamic forces and moments are expressed as linear functions of states. The reduced model is expressed in a form, which only needs the measurable available states. The core objective of the reduced model was to use it in the flight controller design. The reduced model structure appears as follows:

$$\begin{gathered} X=X_0+X_{q}q+X_{\theta }\theta +X_{{\delta }_e}{\delta }_e+X_{{\delta }_f}{\delta }_f,\hspace{1em}Y=Y_0+Y_{P}P+Y_{\varphi }\varphi +Y_{{\delta }_r}{\delta }_r+Y_{{\delta }_f}{\delta }_f, \\ Z=Z_0+Z_{q}q+Z_{\theta }\theta +Z_{{\delta }_e}{\delta }_e+Z_{{\delta }_f}{\delta }_f,\hspace{1em}L=L_0+L_{\theta }\theta +L_{{\delta }_r}{\delta }_r+l_{{\delta }_f}{\delta }_f, \\ M=M_0+M_{\theta }\theta +M_{{\delta }_e}{\delta }_e+M_{{\delta }_f}{\delta }_f,\hspace{1em}N=N_0+N_{\theta }\theta +N_{{\delta }_r}{\delta }_r+N_{{\delta }_f}{\delta }_f. \end{gathered}$$

(4)

The parameter estimation was carried out via the conventional least square estimator. The parameter estimation of the aforesaid forces and moments can be recovered, for example, as follows:

$$\left[\begin{array}{c}{\widehat{X}}_0\\ {\widehat{X}}_q\\ {\widehat{X}}_{\theta }\\ {\widehat{X}}_{{\delta }_e}\\ {\widehat{X}}_{{\delta }_f}\end{array}\right]={\left(R^{T}R\right)}^{-1}{R}^{T}X, R=\left[\begin{array}{ccccc}1& q\left(1\right)& \theta \left(1\right)& {\delta }_e\left(1\right)& {\delta }_f\left(1\right)\\ 1& q\left(2\right)& \theta \left(2\right)& {\delta }_e\left(2\right)& {\delta }_f\left(2\right)\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ 1& q\left(N\right)& \theta \left(N\right)& {\delta }_e\left(N\right)& {\delta }_f\left(N\right)\end{array}\right]$$

(5)

where R is the regression matrix with N observations. The input profile on the elevator to conduct five longitudinal flight maneuvers in the near hover and slow forward flight scenarios included step, doublet, and triplet commands. The parameters were estimated for different flight scenarios. The proximity in the estimated parameter was helpful in choosing a single model for the control design purpose. The estimated parameters are displayed in

Table 1. Average and standard deviation of the estimated parameters for the full and reduced order models [93].

Model	Force	${\mathit{C}}_{\mathbf{0}}$		${\mathit{C}}_{{\mathit{\varphi }}_{\mathit{r}}}$		${\mathit{C}}_{{\mathit{\theta }}_{\mathit{p}}}$		${\mathit{C}}_{{\mathit{\psi }}_{\mathit{y}}}$		${\mathit{C}}_{\mathit{u}}$		${\mathit{C}}_{\mathit{v}}$		${\mathit{C}}_{\mathit{w}}$		${\mathit{C}}_{\mathit{p}}$
	Mom-ent	Avg	Std. Dev.	Avg	Std. Dev.	Avg	Std. Dev.	Avg	Std. Dev.	Avg	Std. Dev.	Avg	Std. Dev.	Avg	Std. Dev.	Avg	Std. Dev.
Full	X	0.1538	0.1131	0.0110	0.0194	0.0483	0.0431	−0.0100	0.0221	−0.049	0.0246	−0.1095	0.1114	−0.0027	0.0358	5.4100E-4	0.0018
	Y	0.0225	0.028	0.0055	0.0316	−0.0205	0.0245	−0.0274	0.0006	8.9522E-4	0.0169	−0.0614	0.0540	−0.0158	0.0331	−0.0073	0.0029
	Z	−0.2861	0.0297	−0.0097	0.0132	0.1699	0.0168	−0.0030	0.0094	0.0014	0.0084	−0.0134	0.0478	−0.0047	0.0121	−4.4929E-4	0.0007
	L	9.6231E-4	0.0004	−8.426E-5	2.970E-5	−6.9285E-4	0.0002	1.4884E-5	3.660E-5	−1.8120E-4	0.0006	7.9049E-5	0.0004	−6.1030E-4	0.0007	−2.1300E-5	5E-6
	M	0.0414	0.0191	0.0033	0.0056	−0.0203	0.0078	0.0041	0.0019	−0.0054	0.0044	0.0014	0.0084	−0.0144	0.0059	−0.0021	0.0006
	N	−0.088	0.0134	−0.0017	0.0012	0.0011	0.0047	−4.3412E-4	0.0011	−2.2278E-4	0.01678	−0.0155	0.0166	0.0011	0.0269	1.2782E-4	0.0002
		${\mathit{C}}_{\mathit{q}}$		${\mathit{C}}_{\mathit{r}}$		${\mathit{C}}_{\mathit{\alpha }}$		${\mathit{C}}_{\mathit{\beta }}$		${\mathit{C}}_{\mathit{\delta }\mathit{f}}$		${\mathit{C}}_{\mathit{\delta }\mathit{a}}$		${\mathit{C}}_{\mathit{\delta }\mathit{r}}$
	X	0.0067	0.0024	0.0067	0.0093	0.0108	0.0019	−0.1252	0.0430	0.0829	0.0037	−0.0035	0.0160	0.0039	0.1243
	Y	3.6187E-4	0.0014	−0.0014	0.0091	0.0265	0.0244	0.0124	0.0237	0.0018	0.0008	−0.0011	0.0021	0.021	0.0008
	Z	−0.0018	0.0010	−0.0051	0.0049	−0.0065	0.0005	0.0519	0.0240	−0.0100	0.0013	4.2333E-4	0.0081	−0.0012	0.0343
	L	7.4992E-5	5.16E-5	1.358E-5	7.86E-6	1.9198E-4	0.0007	1.0647E-4	0.0004	−7.0396E-5	1.18E-6	6.5270E-6	8.14E-5	−5.5306E-6	2.118E-5
	M	0.0013	0.0005	5.516E-4	0.0014	−6.6045E-4	0.0002	−0.0478	0.0181	5.1973E-4	0.0007	−3.1036E-4	0.0035	−1.6626E-4	0.0148
	N	2.7412E-4	0.0008	−0.0019	0.0003	1.9530E-4	0.0245	0.0280	0.0162	0.0033	0.0003	4.2332E-4	0.0032	−2.4090E-4	0.0011
		${\mathit{C}}_{\mathbf{0}}$		${\mathit{C}}_{{\mathit{\varphi }}_{\mathit{r}}}$		${\mathit{C}}_{{\mathit{\theta }}_{\mathit{p}}}$		${\mathit{C}}_{{\mathit{\psi }}_{\mathit{y}}}$		${\mathit{C}}_{\mathit{u}}$		${\mathit{C}}_{\mathit{v}}$		${\mathit{C}}_{\mathit{w}}$		${\mathit{C}}_{\mathit{p}}$
Reduced	X	0.1385	0.0254			0.0354	0.0007
	Y	−0.0353	0.0184					−0.0412	0.0231							−0.007	0.0019
	Z	−0.2437	0.0224			0.1650	0.0022
	L	1.0851E-5	0.0001			9.1882E-6	3.66E-5
	M	−0.0067	0.0191			−0.0031	0.0032
	N	0.0016	0.0057			−7.1896E-4	0.0016
		${\mathit{C}}_{\mathit{q}}$		${\mathit{C}}_{\mathit{r}}$		${\mathit{C}}_{\mathit{\alpha }}$		${\mathit{C}}_{\mathit{\beta }}$.		${\mathit{C}}_{\mathit{\delta }\mathit{f}}$		${\mathit{C}}_{\mathit{\delta }\mathit{a}}$		${\mathit{C}}_{\mathit{\delta }\mathit{r}}$
		0.055	0.0009							−0.017	0.0019	0.0093	0.0023
	Y									0.0028	0.0015			−0.0059	0.0077
	Z	0.0016	0.0014							−0.0012	0.0019	−0.0055	0.0024
	L									−4.0721E-6	5.26E-6			−5.9233E-4	6.68E-6
	M									8.007E-4	0.018	0.0023	0.0039
	N									2.7054E-5	0.0003			0.0208	0.0004

[93].

The model validation was also conducted via two different strategies. In the first method, the validation of the maneuvers was done. In this case, the estimated parameters in the identification cycle were utilized to forecast the aerodynamic forces and moments of other maneuvers. On the other hand, in the second scenario, the validation was done via simulations. In this case, the estimated parameters were used in the nonlinear dynamic simulator, then the states and flight path were reconstructed and compared with the actual flight test data.

The system identification and model validation were carried out via the Grey box method and the parameters were estimated for two different dynamic models, which were named as full order and reduced order for different flight regimes. It was concluded that the developed models were successfully obtained. However, further study is still required to achieve the prediction of the aerodynamic forces and moments with considerable accuracy for a longer period of time and full autonomous flight missions. In this regard, very appealing developments have been reported in [132,133].

3.1.2. System Identification of DelFly II Based on Extended Kalman Filter

Another system identification approach that has been applied to obtain the flight dynamics equations of DelFly II is the extended Kalman filter (EKF) [134]. This approach is the nonlinear version of the Kalman filter that is used to fuse the flight data obtained from the on-board inertial measurement unit (IMU) and external optical tracker OptiTrack [134]. The idea is to use the slow time-scale, but highly accurate data (with the typical average tracking error of 1–2 mm) from the OptiTrack as a basis, especially on cycle-averaged and long-term level, and the fast time-scale data (in this case 512 Hz was used) such as the attitude and the rates and acceleration data from IMU to provide accurate high-frequency data and compensate for error in the tracking data. Therefore, this system identification approach is essentially a closed-loop system identification in which the sensor dynamics are included in the model.

The EKF approach is based on the Kalman filter idea and results in a linear model. However, it is able to deal with nonlinear systems, and it provides optimal estimates to identified parameters. It consists of two phases: the prediction and innovation phases [134]. The main equations are described in the state-space Equations (6)–(12):

$$x(t_{k+1})=f\left(x_k, u_{k, },t_k\right)+g\left(u_k,t_k\right)+w(t_k),$$

(6)

$$z\left(t_k\right)=h\left(x,u,t_k\right)+v\left(t_k\right),$$

(7)

where x, u, z, w, and v denote the states, input, measured output, process noise, and measurement noise, respectively. Both process and measurement noise are assumed to be white and Gaussian, and defined by the covariance matrices Q and R, respectively.

Prediction phase:

$${\widehat{x}}_{k+1,k}={\widehat{x}}_{k,k}+{{\displaystyle \int }}_{t_k}^{t_{k+1}}f\left({\widehat{x}}_{k,k},u_k,t\right)dt$$

(8)

$$P_{k+1,k}={\mathsf{\Phi }}_{k+1,k}{P}_{k,k}{\mathsf{\Phi }}_{k+1,k}^T+{\mathsf{\Gamma }}_{k+1,k}{Q}_{k+1}{\mathsf{\Gamma }}_{k+1,k}^T$$

(9)

Innovation phase:

$$K_{k+1}=P_{k+1,k}{H}_k^T{\left[H_k{P}_{k+1,k}{H}_k^T+R_{k+1}\right]}^{-1},$$

(10)

$${\widehat{x}}_{k+1,k+1}={\widehat{x}}_{k+1,k}+K_{k+1}\left[z_{k+1}-h\left({\widehat{x}}_{k+1,k},u_{k+1}\right)\right],$$

(11)

$$P_{k+1,k+1}=P_{k+1,k}-K_{k+1}{H}_k{P}_{k+1,k}^T$$

(12)

where x_k is the discretized state vector of the linearized system at step k; Φ is the discretized state transition; P is the estimated measurement error covariance matrix; Γ is the corresponding input matrix; K is the Kalman gain; and H is the Jacobian of the measurement equation. Q and R are the covariance matrices that describe the process and measurement noise, respectively.

The EKF was designed to estimate the velocities, body attitude, gyro, and accelerometer biases. The IMU and OptiTrack data were used as the input variables and output measurements, respectively. The OptiTrack provided positions and attitude data, which were computed from the OptiTrack-measured quarternions. From these data, through numerical differentiation and relevant coordinate transforms, the velocities were computed and were fed to the EKF as measurements, $[\begin{array}{ccc}{u}_b^{*}& v_b^{*}& w_b^{*}]\end{array}$. With the z-axis positive upward, the body-frame velocities were given by [134].

$$\left[\begin{array}{c}{u}_b^{*}\\ v_b^{*}\\ w_b^{*}\end{array}\right]=\left[\begin{array}{ccc}\cos \psi \cos \theta & \cos \theta \sin \psi & -sin\theta \\ \cos \psi \sin \varphi -\cos \varphi \sin \psi & \cos \varphi \cos \psi +\sin \varphi \sin \psi \sin \theta & \cos \theta \sin \varphi \\ \sin \varphi \sin \psi +\cos \varphi \cos \psi \sin \theta & \cos \varphi \sin \psi \sin \theta -\cos \psi \sin \varphi & \cos \varphi \cos \theta \end{array}\right]{\left[\begin{array}{c}\dot{x}\\ \dot{y}\\ \dot{z}\end{array}\right]}_{OT}$$

(13)

The quaternions were then used to calculate the Euler angles as it is done in the conventional mathematical modeling of aircraft. The process and measurement equations are given as follows:

Process equations:

$$\begin{array}{l}\dot{\mathsf{\Phi }}=(p-b_p)+(q-b_q)\sin \mathsf{\Phi }\tan \mathsf{\Theta }+(r-b_r)\cos \mathsf{\Phi }\tan \mathsf{\Theta }\\ \dot{\mathsf{\Theta }}=(q-b_q)\cos \mathsf{\Phi }-(r-b_r)\sin \mathsf{\Phi }\\ \dot{\mathsf{\Psi }}=(q-b_q)\sin \mathsf{\Phi }\sec \mathsf{\Theta }+(r-b_r)\cos \mathsf{\Phi }\sec \mathsf{\Theta }\\ {\dot{u}}_B=(r-b_r)v_b-(q-b_q)w_b-g\sin \mathsf{\Theta }+a_z-b_{ax}\\ {\dot{v}}_B=-(r-b_r)u_b+(p-b_p)w_b+g\sin \mathsf{\Phi }\cos \mathsf{\Theta }+a_y-b_{ay}\\ {\dot{w}}_B=(q-b_q)u_b-(p-b_p)v_b+g\cos \mathsf{\Phi }\cos \mathsf{\Theta }+a_z-b_{az}\\ {\dot{b}}_P=0\\ {\dot{b}}_q=0\\ {\dot{b}}_r=0\\ {\dot{b}}_{ax}=0\\ {\dot{b}}_{ay}=0\\ {\dot{b}}_{az}=0\end{array}$$

Measurement equations:

$$\begin{array}{l}{\mathsf{\Phi }}_m=\mathsf{\Phi }+{\upsilon }_{\mathsf{\Phi }}\\ {\mathsf{\Theta }}_m=\mathsf{\Theta }+{\upsilon }_{\theta }\\ {\mathsf{\Psi }}_m=\psi +{\upsilon }_{\mathsf{\Psi }}\\ u_b^{*}=u_b+{\upsilon }_x\\ v_b^{*}=v_b+{\upsilon }_y\\ w_b^{*}=w_b+{\upsilon }_z,\end{array}$$

where x, u, z, v, and w are the states, input, output, process, and measurement noise terms of the filter, respectively, and are defined as follows.

$$\begin{array}{l}x={\left[\begin{array}{cccccccccccc}\mathsf{\Phi }& \mathsf{\Theta }& \mathsf{\Psi }& u_b& v_b& w_b& b_p& b_q& b_r& b_{az}& b_{ay}& b_{az}\end{array}\right]}^T,\\ u={\left[\begin{array}{cccccc}p& q& r& a_z& a_y& a_x\end{array}\right]}^T,\\ z={\left[\begin{array}{cccccc}{\mathsf{\Phi }}_m& {\mathsf{\Theta }}_m& {\mathsf{\Psi }}_m& u_b^{*}& v_b^{*}& w_b^{*}\end{array}\right]}^T,\\ \upsilon ={\left[\begin{array}{cccccc}{\upsilon }_{\mathsf{\Phi }}& {\upsilon }_{\mathsf{\Theta }}& {\upsilon }_{\mathsf{\Psi }}& {\upsilon }_{ub}& {\upsilon }_{vb}& {\upsilon }_{wb}\end{array}\right]}^T,\\ w={\left[\begin{array}{cccccc}{w}_p& w_q& w_r& w_{ax}& w_{ay}& w_{az}\end{array}\right]}^T.\end{array}$$

The states were the attitudes (i.e., the Euler angles $\left[\mathsf{\Phi } \mathsf{\Theta } \mathsf{\Psi }\right])$, body-axis velocities $\left[u_b,v_b,w_b\right],$ and the sensor biases of the gyroscope and accelerometer $\left[b_p,b_q,b_r,b_{ax},b_{ay},b_{az}\right]$. The input was the roll, pitch, and yaw rates $\left[p,q,r\right]$ and the accelerometer-measured accelerations $\left[a_x, a_y, a_z\right]$. The output z were the attitude obtained from the OptiTrack system and velocities. The process noise term v consisted of the measurement noise in the measured velocities and attitudes. Finally, the measurement noise term w consisted of the process noise in the rates and accelerations.

The matrix Q and R were tuned to give more weights to the OptiTrack data except at high frequencies, with Q entries selected to be twice that than was inferred from the data to reduce the confidence in the IMU, and R was 20% smaller than the noise characteristics of the OptiTrack. This choice led to increasingly small innovation errors and faster filter convergence in the EKF. The final covariance matrices are [134]:

$$\begin{gathered} Q=diag\left(3.06, 7.30, 7.12, 44.21, 14.80, 24.20\right), \\ R=diag\left(9.94\times {\mathrm{10}}^{-4},1.07\times {\mathrm{10}}^{-3},8.33\times {\mathrm{10}}^{-4},9.80\times {\mathrm{10}}^{-3},2.47\times {\mathrm{10}}^{-3},6.82\times {\mathrm{10}}^{-3}\right). \end{gathered}$$

Although EKF is well known as a powerful system identification approach that can solve nonlinear estimation problems, it still uses linearized state and/or measurement equations followed by the application of the Kalman filter, and results in a linear estimation problem. Due to the linearization, it yields approximation errors, which tend to underestimate state uncertainties. Therefore, this approach may not be able to capture the nonlinearities and time-varying nature of the dynamics of FWMAVs accurately. This problem can be improved by selecting process equations, which give better descriptions of the dynamics, and in this case, a priori knowledge and engineering insights will be required.

3.1.3. System Identification of Kinkade Slow Hawk Ornithopter

Grauer, in [135], presents the flight dynamics modeling of an avian-based onirthopter that was performed based on nonlinear multi body kinematics and dynamics, which contain the tail and wing aerodynamics and gravitational effects. In this review, we were interested in the system identification approaches, which were performed to obtain the tail and wing aerodynamics through wind tunnel data and free flight data, respectively. The system identification performed was based on a linear least-squares estimation to model the longitudinal aerodynamic forces and torques of the ornithopter. Step-wise regression and equation error in the time domain were utilized to determine the model structure. Once the model structure was determined, parameter estimation was performed using the maximum likelihood cost function,

$$J(\varphi )=\frac{1}{2}{\displaystyle \sum _{i=1}^n{e}^T{W}^{-1}e},$$

where ϕ is the vector that contains all the unknown parameters; e is the residual vector between the outputs and measurements; and W is the inverse noise covariance matrix. This estimator was simplified by the use of equation-error and output-error methods. The equation-error requires process noise assumption, and uses the state derivatives as the measurements, which resulted in a linear estimation problem. Meanwhile, output-error utilized measurement noise and the states as the measurements. This method resulted in a nonlinear and iterative solution that was relatively more accurate than that of the equation-error.

In obtaining the tail aerodynamics, wings were removed and the ornithopter was mounted on a fixed platform. The tail aerodynamics were assumed independent of those of the wings. The system identification was performed using wind tunnel data from four different tests with different free stream velocities and varied angle of attack and sideslip angle. The outcomes of the identification were aerodynamics parameters with a nonlinear model structure. The parameters were obtained through parameter estimation using time-domain equation error [135]:

$$\begin{array}{l}{C}_x=C_{x0}+C_{x_{a^2}}{\alpha }^2,\\ C_y=C_{y\alpha \beta }\alpha \beta ,\\ C_z=C_{z0}+C_{z\alpha }\alpha +C_{z\beta }\beta ,\\ C_l=C_{l_{\beta }}\beta +C_{l_{\alpha \beta }}\alpha \beta ,\\ C_m=C_{m0}+C_{m\alpha }\alpha +C_{m\beta }\beta +C_{m_{{\beta }^2}}{\beta }^2+C_{m_{\alpha V}}\alpha V,\\ C_n=C_{n\beta }\beta .\end{array}$$

The wing aerodynamics were modeled through system identification using flight data, which were obtained using a visual tracking system. The ornithopter, fitted with markers on its body, was flown straight and at a level flight with minimized atmospheric disturbances. System identification was conducted to obtain the mathematical models for the heave force, longitudinal force, and pitching moments that are generated by wings. The generalized aerodynamic forces due to each wing were given as [135],

$$a_{wing}(p,v)={\tau }-[M(p)v+C(p,v)v+g(p)+a_{tail}(p,v)],$$

where p and v are the generalized position and velocity states, respectively; M(p) is a generalized mass matrix; C(p,v) contains nonlinear coupling forces arising from centripetal and Coriolis accelerations; g(p) describes gravitational effects; and a(p,v) describes aerodynamic effects. The forcing vector τ corresponds to forces on the fuselage center of mass, torques on the fuselage center of mass, and torques on the wings. Time-domain equation error and step-wise regression were also utilized in determining the model structure, which was nonlinear in terms of the state variables [135].

$$\begin{array}{l}{C}_x=C_{x0}+C_{x_{{\delta }_w}}{\delta }_w+C_{x_q}q+C_{x{\dot{\delta }}_w}{\dot{\delta }}_w+C_{x_{{\delta }_w^2}}{\delta }_w^2\\ C_z=C_{z0}+C_{z_q}q+C_{z{\dot{\delta }}_w}{\dot{\delta }}_w\\ C_m=C_{m0}+C_{m{\dot{\delta }}_w}{\dot{\delta }}_w\end{array}$$

To handle state variables that were periodic with harmonics of the flapping frequency, a collinearity analysis was performed on the model regressors and followed by an eigensystem analysis. The validation was performed by comparing the aerodynamic forces and moments in the second set of flight data with the predicted ones obtained using the identified models. The heave forces and pitching moment were well predicted, however, the longitudinal force model requires improvement by using longer data records with more airspeed excitation. The discrepancies may be due to the fact that the captured flight data were in fact, closed loop data, which would include the controller and actuator dynamics, and the coupling between the wing and tail aerodynamics that were not well-captured by the model structure.

3.2. Tailless Flapping Wing Micro Air Vehicles

The tailless class of FWMAVs is very appealing compared to the tailed class of FWMAVs because the former provides high agility although it is very sensitive to perturbation and environmental disturbances. Modeling and control of such MAVs is particularly challenging. Therefore, the tailless class has been focused on by a wide number of researchers, and different practical models have appeared [5,10,11,12,14] in the current decade. In the following subsections, we will direct our efforts to highlighting the significant works related to the system identification of tailless FWMAVs.

3.2.1. System Identification of In-Flight Dynamics of Robotic Bee

In this section, the top-down approach to system identification of Robotic Bee [5] is presented. A physics-based model was developed and describes the effects of aerodynamic drag to the translational and rotational dynamics of vertically taking-off FWMAVs. The system identification strategies outlined in [55,126,127,128,136] aimed to quantify the aerodynamic effects and their response to the actuators’ commands in stable air vehicles. However, in unstable systems, closed-loop treatment is carried out to keep away the undesirable coupling between input commands and dynamic responses in [137,138]. In the context of close loop identification, the Grey box model identification is one of the suitable strategies [138]. In order to carry out system identification, the original nonlinear model is presented in such a way that the vector of unknown parameters can be evaluated via linear regression. This strategy provides certain benefits such as global extremum and analytic solutions. This helps in eliminating the tuning of the optimizing parameters and eliminating the need of a local search that may increase the relative computational complexity [139]. At this point, the in-flight aerodynamic model could be obtained by approximating the instantaneous drag, D, on the two wings, with high Reynolds number law [140,141].

$$D=2\times \frac{1}{2}\rho S{C}_d{v}^2,$$

(14)

where ρ is the density of the fluid; S is the surface area of the wing; $C_d$ is the coefficient of drag; and v is the speed of the wing relative to air. In the case of deviation of the robot from hovering and with additional damping force, the drag force contributes to the rotational damping, thus, the alternate form of (14) is [5]

$$D\approx mbu+m\alpha \omega$$

(15)

where m is the mass of the robot; b is the normalized drag force; u is the speed in forward direction; $\alpha$ is the angle of attack; and ω is the rotational speed of the robot.

The resulting drag moment can be expressed as follows:

$$M\approx Jc\omega +J\beta u,$$

(16)

where c and β are the rotational damping coefficients normalized by inertia. Due to (15) and (16), the viscous forces and torque from the wings, which dominate the drag forces and are shown critical to the open loop dynamics and stability of FWMAVs [34,141,142], can also be treated as linear functions.

The translational dynamic equation of the robot can be modeled as follows:

$$m\ddot{X}=m\mathsf{\Gamma }\widehat{Z}+W+F_a,$$

(17)

where $F_a$ represents the lumped aerodynamic forces; $W=mg$ represents the gravitational force; Γ is the normalized thrust with dimension of acceleration; and $\widehat{Z}$ represents the axis along the body of the robots. It includes the additional aerodynamic force due to drag forces on the wings and viscous forces on the body that can be further modeled as linear functions of translational and angular velocities. Equation (17), with explicit expression of $F_a$, can be re-expressed as follows:

$$\frac{d}{dt}\dot{X}+\omega \times \dot{X}=\left[\begin{array}{c}0\\ 0\\ \mathsf{\Gamma }\end{array}\right]-g\left[\begin{array}{c}{R}_{31}\\ R_{32}\\ R_{32}\end{array}\right]+\left[\begin{array}{c}-b_x\dot{x}-{\alpha }_x{\omega }_y\\ -b_y\dot{y}+{\alpha }_y{\omega }_x\\ -b_z\dot{z}\end{array}\right],$$

(18)

where b_i and α_i are contributed from the body velocity and angular velocity with i ∈ {x, y, z}. When the FWMAV is not in perfect hovering, then it is affected by additional torques, therefore, the resultant torque on the robot can be expressed by the following rotational dynamic equations:

$${\displaystyle \sum }\tau =J\frac{d\omega }{dt}+\omega \times J\omega ={\tau }_c+{\tau }_a+J_d,$$

(19)

where $J_d$ is the unknown normalized disturbance. Note that in this presentation, the roll and pitch were focused and the aerodynamic contributions were modeled as linear in its variables. In more detail, this equation can be expressed as follows:

$$\frac{d}{dt}\left[\begin{array}{c}{\omega }_x\\ {\omega }_y\end{array}\right]=\left[\begin{array}{c}{\tau }_{c, x}{J}_x^{-1}\\ {\tau }_{c, y}{J}_y^{-1}\end{array}\right]+\left[\begin{array}{cc}\left(J_y-J_z\right)J_x^{-1}& 0\\ 0& \left(J_z-J_x\right)J_y^{-1}\end{array}\right]\left[\begin{array}{c}{\omega }_x{\omega }_z\\ {\omega }_y{\omega }_z\end{array}\right]+\left[\begin{array}{c}-c_x{\omega }_x+{\beta }_x\dot{y}\\ -c_y{\omega }_y-{\beta }_y\dot{x}\end{array}\right]+\left[\begin{array}{c}{d}_x\\ d_y\end{array}\right].$$

(20)

The aforementioned Equations (18) and (20) are linear in their structures and therefore, the estimation of unknown parameters via the linear regression approach is the ultimate objective. The linear regression model could be expressed as follows:

$${\mathsf{\Psi }}_i={\phi }_0+{{\displaystyle \sum }}_{j=1}^n{\chi }_{ij}{\phi }_j+{ϵ}_i; i\in \{1,2,\dots ,N\},$$

(21)

where ${\phi }_0$ represents an unknown possible offset or affine term; ${ϵ}_i$ represents the residual of an $i^{th}$ observation; N is the number of total observations; n is the number of unknown parameters; and Ψ is an N × 1 vector, which needs to be extracted. The following least square estimator is considered to provide the estimated parameters.

$$\widehat{\phi }={({\chi }^T\chi )}^{-1}{\chi }^T\mathsf{\Psi }, \mathrm{and} R^2=1-\frac{{\sum }_{i=1}^N{ϵ}_i{}^2}{{\sum }_{i=1}^N{\left({\mathsf{\Psi }}_i-\frac{1}{N}{\sum }_{j=1}^N{\mathsf{\Psi }}_j\right)}^2}$$

(22)

where $R^2$ indicates the goodness of the estimation using curve fitting with $R^2$ close to 1 for a good estimation.

The translational dynamics presented in (18) in the body fixed frame could be further subdivided into lateral and altitude dynamics. The estimates of the lateral dynamics into a merged equation can be expressed as

$${\mathsf{\Psi }}_i=\{\begin{array}{cc}\ddot{x}+{\omega }_y\dot{z}-{\omega }_z\dot{y}& {\mathrm{for} \mathsf{\chi }}_i=\left[\begin{array}{ccccc}-R_{31}& -\dot{x}& -{\omega }_y& 0& 0\end{array}\right] \\ \ddot{y}+{\omega }_z\dot{x}-{\omega }_x\dot{z} & {\mathrm{for} \mathsf{\chi }}_i=\left[\begin{array}{ccccc}-R_{32}& 0& 0& -\dot{y}& -{\omega }_x\end{array}\right]\end{array}$$

(23)

and $\phi ={\left[\begin{array}{ccccc}g& b_x& {\alpha }_x& b_y& {\alpha }_y\end{array}\right]}^T$ as the vector of parameters to be estimated. Furthermore, for altitude dynamics, it is expressed as follows

$${\mathsf{\Psi }}_i=\{\begin{array}{cc}\ddot{z}+{\omega }_x\dot{y}-{\omega }_y\dot{x}& , \mathrm{for} {\mathsf{\chi }}_i=\left[\begin{array}{ccc}\frac{T}{s{\gamma }^{-1}+1}& -R_{33}& -\dot{y}\end{array}\right]\end{array}$$

(24)

and $\phi ={\left[\begin{array}{ccc}\mu & g& b_z\end{array}\right]}^T$ as the vector of parameters to be determined. T is the commanded thrust and s is a Laplace variable.

The rotational dynamics, because of the existence of coupling of the moment of inertia between the axes, are less sophisticated. These terms are lumped into ${\nu }_i$ and $k_i$ variables. The term ${\nu }_i$, for example, indicates the inverse of the moments multiplied by some scaling factor, which represents the relation between the command torque and the actual torque. The factor ${\left(s{\gamma }^{-1}+1\right)}^{-1}$ indicates the response time of the system, and is very much similar to the response time of the system to the thrust command. One may have the following expressions

$${\mathsf{\Psi }}_i=\{\begin{array}{cc}\frac{d}{dt}{\omega }_x ,& { \mathsf{\chi }}_i=\left[\begin{array}{cccc}\frac{{\tau }_{c, x}}{s{\gamma }^{-1}+1}& {\omega }_y{\omega }_z& -{\omega }_x& \dot{y}\end{array}\right]\end{array}\},$$

(25)

and $\phi ={\left[\begin{array}{cccc}{\nu }_x& k_x& c_x& {\beta }_x\end{array}\right]}^T$ as the vector of unknown parameters, and

$${\mathsf{\Psi }}_i=\{\begin{array}{cc}\frac{d}{dt}{\omega }_y& {\mathsf{\chi }}_i=\left[\begin{array}{cccc}\frac{{\tau }_{c, y}}{s{\gamma }^{-1}+1}& -{\omega }_x{\omega }_z& -{\omega }_y& -\dot{x}\end{array}\right]\end{array}\},$$

(26)

and $\phi ={\left[\begin{array}{cccc}{\nu }_y& k_y& c_y& {\beta }_y\end{array}\right]}^T$. The available data from the flight maneuvers were downsampled using cubic spline to improve the continuity and a 3rd order Chebyshev-type-II low pass filter was used to smooth the trajectories. The cutoff frequency of the filter was 50 Hz. The angular velocities were calculated via the use of rotation matrices and missing states were reconstructed via differentiation. The body fixed frame data were recovered via projecting the available states into the body coordinate frame. The estimated parameters were categorized into three for lateral dynamics, altitude dynamics, and rotational dynamics.

3.2.2. System Identification via Frequency Approach of Drosphilla melanogaster

Frequency-based system identification methodologies were studied in [116,117,142,143]. In this section, the methodology in [143] for the estimation of stability and control derivatives is outlined. Since the flight dynamic equations, in general, are very complex, the original quasi-steady model proposes the lift and drag as functions of the instantaneous wing motion. Therefore, some basic assumptions are made to proceed to the controllability analysis. The first assumption portrays the quasi-steady model [88] as an experimental fit to the scaled dynamics of the flapping flier. However, the actual dynamics contain the lift and drag forces as functions of instantaneous wing motion, the aerodynamic forces are functions of the body motion, and the mechanism was the focus of the study in [142,144,145,146], which gave birth to passive aerodynamic damping. The second assumption is that the instantaneous lift and drag forces on the rigid body are approximated by wing stroke average forces, and furnish the stability and control derivatives of the aircraft. In this scenario, the time dependency is removed via the averaging theory [147]. Consequently, a linear time-invariant (LTI) system of the following form is equipped.

$$\dot{x}=Ax+Bu,$$

(27)

where $A\in R^{n\times n}$ represents the state matrix; $B\in R^{n\times p}$ is the control derivative matrix; and $u\in L_2^p\left[0, \infty \right)$ indicates the controlled input history and $x\in L_2^n\left[0,\infty \right)$ points to the state history of the system. The aforementioned components smoothen the way to system identification. To obtain stability and control derivatives and to perform nonlinear simulations, the system was considered as a rigid body in combination with quasi-steady aerodynamics and the state perturbations caused by the state deviation from the equilibrium. Similar work can be found in [142] where the quasi-steady aerodynamics formulation was treated as state perturbation from a reference equilibrium condition via proper definition of the wing tip speed and the wing angle attack. The current subsection focused on the study in [143], with the aerodynamic model outlined in [142], which included the effects of the rotational and translational rigid body motions.

In the underlying strategy, the trim kinematics [99,142] were treated as periodic, represented by ϕ(t), in a stroke plane, which makes an angle $\beta$ with the horizontal. Both the left and right wing were actuated simultaneously with the same control input for the said longitudinal motions. The control inputs were defined as stroke plane inclination ${\beta }_c=\frac{1}{2}\left({\beta }_R+{\beta }_L\right)$ [117,148], stroke plane offset ${\varphi }_{0ff}=\frac{1}{2}\left({\varphi }_{0f{f}_R}+{\varphi }_{0f{f}_L}\right)$ [143,148], and asymmetric wing angle ${\alpha }_{ad}=\frac{1}{2}\left({\alpha }_u-{\alpha }_d\right)$ [117,143,148]. These inputs are shown in Figure 5.

Note that the input ${\beta }_c$, which physically represents the tilting of the stroke plane, is used to generate the pitch moment and forward force [142], ${\varphi }_{off}$ physically represents the force shift of the wing sweep for pitch moment generation and ${\alpha }_{au}$ in a similar way represents the upstroke or downstroke asymmetry in wing angle relative to the stroke plane. This input is used to generate the forward force.

The longitudinal dynamics of the insect were excited by applying the chirp signal at each input. The overall transfer function was found via the spectral analysis of the input transfer function $G_I\left(\omega \right)$ and the output transfer function $G_O\left(\omega \right)$.

$$G_s\left(\omega \right)=\frac{G_O\left(\omega \right)}{G_I\left(\omega \right)},$$

(28)

For the sake of completion, the input/output response in terms of magnitude, phase, and coherence for ${\varphi }_{off}$ to θ is shown in Figure 6. The haltere-based feedback on the pitch axis was used to make the system matrix Hurwitz. The system identification strategy via the frequency-based analysis is considered to be more powerful [142] than the analytical linearization and computational based techniques because by coherence (γ(ω)), one may easily determine the model over certain frequencies where the linear and nonlinear model closely match.

$$\gamma \left(\omega \right)=\frac{{\left|G_{IO}\right|}^2}{G_I\left(\omega \right)G_O\left(\omega \right)}$$

(29)

The stability and control derivatives, in the forthcoming general state space model, were determined at low frequency regions (0–20 Hz) by maximizing the coherence. In addition, the small periodic high frequency motions were neglected in this range of frequencies. The resulting mathematical model for longitudinal dynamics is shown by Equation (30).

$$\left[\begin{array}{c}\mathsf{\Delta }\dot{u}\\ \mathsf{\Delta }\dot{w}\\ \mathsf{\Delta }\dot{q}\\ \mathsf{\Delta }\dot{\theta }\end{array}\right]=\left[\begin{array}{cccc}{X}_u& 0& 0& -g\\ 0& X_w& 0& 0\\ M_u& 0& M_q& 0\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{c}\mathsf{\Delta }\mathrm{u}\\ \mathsf{\Delta }w\\ \mathsf{\Delta }\mathrm{q}\\ \mathsf{\Delta }\mathsf{\theta }\end{array}\right]+\left[\begin{array}{cccc}0& X_{\beta }& X_{{\varphi }_{0ff}}& X_{{\alpha }_{ad}}\\ Z_{\mathsf{\Phi }}& 0& 0& 0\\ 0& M_u& M_u& M_u\\ 0& 0& 0& 0\end{array}\right]\left[\begin{array}{c}{\mathsf{\Phi }}_c\\ {\beta }_c\\ {\varphi }_{0ff}\\ {\alpha }_{ad}\end{array}\right]$$

(30)

The typical determined values are displayed in

Table 2. This indicates the estimated stability and control derivatives [143].

Parameter	Value	Standard Deviation
Xu	−12.3	4.0%
Zw	−4.7	9.1%
Mu	547	3.6%
$M_q$	−33.3	4.6%
$X_{{\varphi }_{off}}$	12.7	12.2%
$X_{\beta }$	20.7	7.6%
$Z_{\varphi }$	−17.4	4.8%
$M_{{\varphi }_{off}}$	6028	3.6%
$M_{\beta }$	−2826	3.8%
$X_{\alpha }{}_{ud}$	−12360	5.0%
$M_{\alpha }{}_{ud}$	−87.74	10.7%

where the standard deviations were measured by following the procedure of [149]. In the state space model, an additional control input ${\varphi }_c$ was introduced, which physically indicates the stroke amplitude and which has a direct effect on the heave dynamic without affecting the other system states.

The discussed frequency-based study has far appeal with fruitful results. In the subsequent subsection, we highlight the individual efforts for the explicit identification of either the system’s stability derivatives or system’s stability derivatives and system’s control derivatives.

3.2.3. Computational Fluid Dynamic Approach Based System Identification

The aerodynamics and system modeling of the insects and insect mimicking FWMAVs have attracted a wide number of researchers [34,46,84,85,150,151,152,153,154]. In recent developments, the dynamic flight stability of insects has gained much attention [86]. They focused on a formal quantitative study regarding the dynamic flight stability of desert locust Schistocerca Gregaria. In their work, they estimated the stability derivatives via experimental study that were further studied for eigenvalue analysis to recognize different stability modes like the unstable oscillatory mode, stable fast subsidence, and slow subsidence mode. Note that the experiment carried out was subjected to some limitations such as that the insect had to be tethered and the reference flight might not have been in equilibrium. In practice, the presence of such limitations may make the estimation of some morphological parameters impossible. In order to pave the way to estimate stability derivatives for all conditions and different flight maneuvers, [112] presented a CFD approach to measure the stability derivatives of a hovering bumblebee while capturing the flight data used in the study of [84,85], and the stability derivatives of beetle-mimicking FWMAV in [155]. In this study, they first determined the equilibrium condition, while taking account of the rigid body assumption, and utilized the flow equations and solution methodology presented in [48] to compute the stability derivatives for the longitudinal flight dynamics. In the stability analysis, once again, they identified the same mode as were concluded in [86]. In [156], the work of [112] was extended to carry out a formal analysis of the stability, controllability, and the stabilizing control of the hovering insect. The stability analysis resulted in the same modes as identified in [86]. Therefore, a control action in the feedback loop was needed to stabilize the unstable oscillatory mode, and stability augmentation was pursued for the slow subsidence mode. The stabilization was done by feeding back a controlled input (i.e., the changes in the mean stroke angle as a function of pitch attitude, pitch rate, and horizontal velocity), whereas the second case was handled via a feedback control input (i.e., the changes in the stroke amplitude as a function of vertical velocity). This study was fascinating for further study in different flight envelopes. In the context of longitudinal stability analysis, in [111], the dynamic stability of the beetle-mimicking FWMAVs [43] was investigated. The authors followed the trim condition determination strategy of [86], and CFD based stability derivatives were computed. The stability analysis was concluded as very similar to [86,111]. In order to determine the flight stability, one needs to have the identified dynamics for both longitudinal and lateral motions. Currently, a quantitative analysis for the lateral dynamic stability of a hovering insect is addressed in [157] via the CFD strategy [111,156]. In this case, they focused on the drone fly instead of the bumblebee. They have computed the stability derivatives and, once again, the eigenvalue-based analysis was conducted, and very similar results to the longitudinal scenario (i.e., one unstable slow divergence mode, one stable slow oscillatory mode, and a stable fast subsidence mode) were identified. The aforementioned CFD strategies provide a suitable understanding about the stability analyses of FWMAVs. However, the CFD-based approach is computationally costly and time-consuming.

3.2.4. Parameter Estimation for Linear Damping Model of DelFly Nimble

The tailless version of DelFly, known as DelFly Nimble, was used in closed-loop modeling of the longitudinal dynamics of a tailless flapping wing robot [158]. The idea of the modeling is to investigate the effect of changing the robot physical parameters and the controller gains to the robot’s flight dynamics. Due to the lack of tail that provides passive stabilization, the robot relies on the wing actuation system and active control system to provide its stability and control. The linear damping model was selected as it can approximate the aerodynamics of complex flapping wings with a minimal number of parameters in hovering and forward cruise flight (with quasi-steady assumption) that were measured in a wind tunnel, as shown in [139,159]. It has also proven its ability in predicting the horizontal drag forces of flapping wings. In [158], vertical dynamics and pitch control mechanisms were included. The overall model structure consists of controller dynamics and open loop dynamics, which are comprised of the actuator dynamics and the longitudinal flight dynamics.

System identification was used to determine the dihedral actuator dynamics and force coefficients, $b_x$ and $b_z$, which define the longitudinal and vertical components of the aerodynamic drag forces as described by the following equations.

$$\begin{array}{l}{D}_x=-b_{x}fu,\\ D_z=-b_{z}fw,\end{array}$$

where f is the flapping frequency; u and w are the longitudinal and vertical body velocities, respectively. Flight data were collected as the robot was flown in an open jet wind tunnel at various wind speeds ranging from 1 ms⁻¹ to 2.4 ms⁻¹. Both u and w were computed based on the averaged wind speed and body pitch. Linear least squares method was used to determine $b_x$ and $b_z$. Validation was performed by comparing the estimated state derivatives from the model with those derived directly from the flight data. The overall results show that the developed model can capture the robot’s dominant dynamic effects and can be useful for initial tuning of the control parameters. The performance of this approach, however, will depend on the controller used.

4. Conclusions

A comprehensive review on system identification and mathematical modeling of flapping wing micro-aerial vehicles has been presented. The details were provided for tailless and tailed natural fliers and micro-aerial vehicles such as Drosophila, hawk ornithopter, RoboBee, KUBettle, and Delfly-II as well as the Delfly Nimble. The outcomes of the system identification approaches in all reviewed articles were satisfactory in general. The obtained models were validated to predict and/or estimate the dynamics of the plants used, and they could achieve reasonable accuracies. The identification strategies that have been used were top-down Grey box approaches. These include the extended Kalman Filter, linear regression and least square estimation, frequency-based identification, and computational fluid dynamic based system identification. Most of the approaches were used to estimate model parameters that are related to aerodynamics and result in standard aircraft mathematical models or equations of motion, which are adopted with or without modifications. The modifications were based on a priori knowledge of the natural fliers’ aerodynamics and engineering insights. The final models typically involve decoupling of the longitudinal and lateral dynamics. Additionally, wind tunnel or CFD-based system identification has been mostly performed to determine the stability derivatives, which were then again used together with the conventional aircraft equations of motion. CFD-based procedures were performed on the tested MAVs, while they were either mounted on a fixed platform or were tethered, thus, the natural flight motions of the vehicles were constrained, and this would change the aerodynamics during the flight. The final mathematical model would be either a full-model or a reduced model, in which insignificant stability derivatives were removed. Moreover, we found that a more complete system identification process aimed for design decision and controller development was presented in [159], in which a simplified linear damping model of the DelFly Nimble was obtained.

This review has shown that system identification approaches have a promising potential to uncover many unanswered questions in generating valid mathematical models of FWMAVs that can capture the overall FWMAVs’ complex flight dynamics. From the review, we conclude that system identification approaches have not been used to their maximum advantage for modeling FWMAVs at the present stage. The approaches were mostly used to obtain the aerodynamic parameters of certain FWMAVs for hovering and/or steady flight, linear regression and least square estimation were mostly used, and the acquired models were mostly linear. Thus, they are valid only for particular flight modes. There are still questions on the fidelity of the models obtained, particularly those for various maneuvers.

Better outcomes of system identification and modeling complex systems such as FWMAVs can be achieved with detailed analyses done before the start of the identification that lead to appropriate model selection. Such analysis will be critical to determine the optimized model structure and model order where a priori knowledge and engineering insights are combined with formal properties of models. In addition, it is imperative to have clear objectives of the system identification to be performed, for example, system identification for design decision, simulation or validation, and/or for flight control design. Therefore, specific experimental data, efficient model selection, and structure as well as identification approach can be determined accordingly. Many nonlinear system identification approaches that are known for their powerful capabilities in dealing with nonlinear dynamical systems such as the unscented Kalman filter (UKF) and nonlinear autoregressive moving average model with exogenous inputs (NARMAX) have not been considered to solve the modeling problems of FWMAVs. Together with a priori knowledge and engineering insights, those approaches are expected to provide more accurate outcomes than those of the linear ones.

Furthermore, there have hardly been any black-box approaches such as neural-network based system identification reported in use for the system identification of FWMAVs. Such approaches have been used in the identification of other types of unmanned aerial vehicles such as in [160]. In conclusion, more sophisticated system identification strategies, combined with knowledge of aerodynamics and flight dynamics of natural flapping wing fliers, should be utilized, so that we can gain more understanding of highly complex dynamics of FWMAVs, particularly to achieve the autonomous flight of FWMAVs, and develop energy storage mechanisms to obtain energy efficient flight of insect fliers [161].

5. Future Works

Finally, this review article is furnished via our understudy theoretical strategy for system identification. This strategy followed the Grey box approach while utilizing the differentiator observer theory [162]. We proposed a sliding mode-based system identification that utilizes the robust and accurate nature of sliding mode observer and flatness property. The proposed technique is briefly discussed in this section. The motivation for the use of the aforementioned technique is its inherent robustness in the face of uncertain conditions. In this work, we also relax the identifiability condition [163]. At this point, the following assumptions are furnished.

• It is assumed that the control input derivatives are known and certain.
• The system’s stability derivatives are assumed to be unknown and separable.
• The states of the system are observable.
• The system under study must satisfy the flatness property [164].

Before proceeding to the stability derivative estimation, the observability and separability conditions are checked offline. The flatness property plays a significant role in this strategy based on which parameters are extracted. The controller used in the closed-loop estimation is the one designed via sliding mode control strategy [165]. The identification strategy flowchart is shown in Figure 7. At the first step, the mathematical model of the system is developed and some parameters are marked as unknown in an uncertain environment. Then, by utilizing the flatness property, these parameters are extracted and are further estimated by invoking the known data. The mathematical model is finalized for control design purpose by using these estimated parameters. Therefore, we can consider our strategy as the parameter identification of the flight dynamics of the FWMAVs.